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COMPUTER PROGRAM FOR FINITE -DIFFERENCE SOLUTIONS 


OF SHELLS OF REVOLUTION UNDER ASYMMETRIC 
DYNAMIC LOADING 

By Wendell B. Stephens and Martha P. Robinson 
Langley Research Center 

SUMMARY 

A general computer program written in FORTRAN IV language which determines 
the linear asymmetric bending behavior of a statically or dynamically loaded elastic thin 
shell of revolution is presented. The loading may be applied either mechanically or 
thermally. The variables are separated by representing the loads, displacements, and 
stresses by Fourier series expansions in the circumferential direction. The resulting 
set of equations is solved numerically by using finite -difference approximations in the 
meridional direction and backward differences in the time direction. A three-layered 
cross section which is symmetric about the middle surface is allowed. The boundary 
conditions are taken in a general form which allows the program to handle elastic 
restraints specifying a linear combination of edge forces and displacements. Nonhomo- 
geneous initial conditions are also allowed. The data input procedure is described in 
detail and sample calculations are included. 

INTRODUCTION 

The linear elastic behavior of any shell of revolution with a static asymmetric load 
has been programed in reference 1 by using Fourier series expansions along the circum- 
ference along the meridian of the shell. The programed analysis contained in refer- 
ence 1 is based on the analytic formulation presented in reference 2. The shell theory 
used is that of reference 3. The program in reference 1 calculates the Fourier coeffi- 
cients of the series but does not perform the summations of the series. This paper 
extends the analytic formulation and computer program of reference 1 to include asym- 
metric dynamic response of shells of revolution and includes provisions for summation 
of the terms of the series. Numerical integration of the dynamic equations is similar to 
that given in reference 3 for a cylindrical shell and is based on Houbolt's backward dif- 
ference method (refs. 4 and 5). The loading on the shell may be either mechanical or 
thermal and these loads may be either static or dynamic. The thickness of the shell may 
vary along the meridian, but the shell cross section must be symmetric about the shell 


middle surface. In addition, the initial conditions include provisions for initial deforma- 
tions. This paper is a user’s document which contains the necessary instructions for 
preparation of input data and subprograms. The program is written in the Control Data 
version (ref. 6) of FORTRAN IV language for operation in the scope 3.0 digital computer. 
The program requires an octal storage of 70 000 memory words. The output of the pro- 
gram lists the shell description and the nondimensional Fourier series summations of 
the displacements, rotations, moments, and force resultants in a tabular (columnwise) 
format. 

The program presented in this paper is illustrated by two dynamic -response 
example problems. In the first problem a comparison is made with exact results for 
axisymmetric deformation of a cylindrical shell under initial deformation. The second 
example demonstrates the input preparation for a conical shell with asymmetric deforma- 
tions. The second example problem is a planetary entry "aeroshell" subject to asym- 
metric pressures as the shell passes through the atmosphere with a small oscillation. 

The data preparation is discussed in detail for this practical application. 

SYMBOLS 

The units used for the physical quantities in this paper are given both in the U.S. 
Customary Units and in the International System of Units (SI). Factors relating the two 
systems are given in reference 7 and those used in the present investigation are pre- 
sented in appendix A. 

a reference (or characteristic) length 

b nondimensional membrane stiffness defined in appendix C 

d nondimensional bending stiffness defined in appendix C 

E 0 reference modulus of elasticity 

f^ n ) nondimensional transverse meridional shear (see appendix B) 

f frequency 

h shell thickness 

h 0 reference thickness 
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last boundary station 
bending-moment resultants 

nondimensional thermal-moment resultant defined in appendix C 

modified twisting moment (eq. (8)) 

nondimensional Fourier coefficients for bending moments 
(see appendix B) 

total number of stations 

membrane force resultants 

modified membrane shear (eq. (7)) 

Fourier index 

first boundary station 

Fourier coefficients for loads (see appendix B) 
transverse shear resultants 

distributed loads in normal, meridional, and circumferential directions, 
respectively 

aerodynamic pressure 

radial distance from axis of symmetry to shell middle surface 
shell meridian 

Fourier coefficient for temperature 

Fourier coefficient for midplane temperature variation (eq. (30)) 
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Fourier coefficient for temperature gradient per unit thickness normal 
to middle surface (eq. (30)) 

inverse of frequency f; period 

face sheet thickness 

time 

nondimensional Fourier coefficients for membrane force resultants 
nondimensional thermal -force resultant defined in appendix C 
meridional and circumferential displacements 

nondimensional Fourier coefficients for meridional and circumferential 
displacements (see appendix B) 

normal displacement 

nondimensional Fourier coefficient for normal displacement 
(see appendix B) 

coefficient of thermal expansion 

angle between a normal to the shell and fluid flow 

coefficients of backward difference expression 


y = p t /p 


A 


e 


meridional increment between interior stations 


time increment, e = / — AT 

Vpa^ 


coordinate normal to and orginating at midsurface of shell, positive 
outward 


circumferential coordinate 
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Matrices: 


A,B,C,D,E, ^1 
F,G,H,J,fi,Aj 

e, z,y/l 

f, g/ j 


Poisson's ratio 

nondimensional radius, r/a 
shell density 
density of face sheets 
reference stress 
nondimensional time (eq. (9)) 

meridional rotation 
colatitude angle 

nondimensional meridional rotation (see appendix B) 
angle of attack 

amplitude of aeroshell oscillation (eq. (36)) 
nondimensional curvatures (eqs. (1) and (2)) 


4X4 matrices 


4x1 matrices 


A prime indicates a derivative with respect to £. A dot indicates a derivative 
with respect to r. Note that the superscript n is dropped from the Fourier coeffi- 
cients when doing so would not cause confusion. 
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ANALYTICAL FORMULATION 
Shell Geometry 

The shell geometry and coordinates are shown in figure 1 and are identical to 
those used in reference 1. A point on the shell is specified by coordinates where 

£ = s/a is the nondimensional meridional coordinate, s is the meridional coordinate, 
a is a reference dimension of the shell, 6 is the circumferential coordinate, and 
£ is a coordinate normal to and originating at the middle surface, positive outward. 




Figure 1.- Surface geometry and coordinates. 
OP = rj O^P = a/cuj; and O 2 P = a/o>g. 


If the shape of the middle surface is given by p = p(£) where p = r/a and r is the 
distance OP, the nondimensional principal curvatures can be written as 


[i - (p') 2 1 1/2 

P 


(1) 




y' + 




( 2 ) 


6 



where 


and the prime indicates a differentiation with respect to 


( 3 ) 


Dynamic Response Terms 

The equations of motion for a shell are obtained from the equilibrium equations of 
reference 1 by the addition of the acceleration terms as follows: 

+ ^ - p'm 0 J + af^y + 4 N |0 - p N 0 j 
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where and are defined in reference 3 as 


n £ 6 > = |( N 40 + N 04 ) + " W |)( M £0 " M 0 


l) 
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40 -|( m 40 + m 04) 


The nondimensional time r is 


( 8 ) 



..At 
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The positive directions of the forces, moments, rotations, and displacements are shown 
in figure 2. 



(c) Displacements. 



(d) Moment resultants. (e) Loads per unit area. 

Figure 2.- Positive sense of force, moments, shears, displacements, 
and loads on a shell segment. 

The Fourier series expansions of the unknowns in equations (4) to (8) and of <t>| 
and are summarized in appendix B in appropriate nondimensional form. By uti- 

lizing these relationships and the stress-strain and strain-displacement relationships, 
the governing second-order partial differential equations in matrix notation become 

Ez^ n ) + Fz( n ) + Gz( n ) = e + Dz( n ) (10) 

where 
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( 11 ) 



V J 


and u^ n >, U W and are amplitudes of the Fourier harmonics of the 

three displacements and the meridional moment. The superscript n refers to the 
nth Fourier harmonic and is omitted hereafter for purposes of simplification in notation. 
The dots represent derivatives with respect to time r. The elements of the E, F, G, 
and e matrices are defined in both references 1 and 2 and are included in appendix C. 
The D matrix is 


1 



0 

0 

0 


0 

1 

0 

0 


0 0 

0 0 

1 0 

0 °_ 


The vector y of stresses and moments is related to the vector z of displacements 
and rotations by the equation 


y - Hz'+ Jz + f 


( 12 ) 


where the H, J, and f matrices are defined in reference 2 and appendix C. The 
vector y is defined to be 


y = 



v j 


(13) 


The t|, t^ fl , f^, and components are the Fourier coefficients (see appendix B) 
of the axial -stress resultant, transverse -stress resultant, shear stress, and rotation 
variables. It is convenient to express the boundary conditions at either or £ 
the end points of the meridional generatrix, in terms of the equation 
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S2y + Az = l 


(14) 


where the £2 and A matrices and l vector are chosen so that the prescribed dis- 
placements and forces at the boundary are satisfied. It should be noted that the dynamic 
inertia terms affect only the equilibrium equations but not the boundary conditions. 

Computational Techniques 

The central finite-difference approximations used in equation (10) along the 
interior stations of the meridian are defined in reference 1 as 


„ ' _ z i+l z i-l 
1 " 2A 


(15) 


T T 


z i+l ^ z i + zi_i 
A2 


(16) 


where 


A = 


N - 2 


(17) 


and where i = 2, 3, . . ., N- 1. Here N is the number of meridional stations. The 
first derivative of z and value of z at the boundaries of the meridian are 


z 0 = f( z l +z 2) 
z o’ = i(“ z l + z 2) 

Z L = §( Z N-1 +z n) 
Z L = a ( _Z N-1 + Z N 


) 


(18) 


The first shell edge is located midway between stations 1 and 2 and the end of the merid- 
ian is located midway between stations N - 1 and N. The first and final shell edges 
are denoted by subscripts 0 and L, respectively. 
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The time derivatives can be represented by backward differences (as in refs. 5 
and 6): 


h ,1 = ®l z i,l + Pl z i,l-1 + n z i,l-2 + h z i,l-3 + a l z i,0 + h' z i,0 

(f = 0,1,2, ... , i = 1,2, . . .,N) (19) 

where dTj , /% , y z , and 6 \ are constants which depend on the time step and which, in 
general (l S3), define a four-point backward second derivative. The constants and 

/3 1 are required for including nonhomogeneous initial conditions z^ q and z^q. The 
first subscript on z denotes the spatial station and the second subscript denotes the 
time station. Since z^ q and z^ q are given initial conditions that allow z^ q to be 
calculated from equation (10), fictitious time points at l = - 1 and l = - 2 can be 
obtained by using the difference expression 


z i,0 = ^( 2z i,l + 3z i,0 - 6z i,-l + z i, -2) 


( 20 ) 


z i>0 = To ( z i,l ” 2z i,0 + z i,-l) 


( 21 ) 


where e is a time increment of r. Thus equations (20) and (21), the given initial condi- 
tions z^ q and z^ q, and equations (10) and (12) in finite -difference form comprise 
enough information to define the coefficients of equation (19). Therefore, at 1=0, 


<*0 = 00 = >"0 = 6 0 = “0 = 0 

^0 = 1 


( 22 ) 


at 1=1, 






"\ 


yj = 6j = 0 ) 
6 


“1=F 


/3i = -2 


(23) 
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at 1=2, 



Equation (19) with the constants in equation (25) is the standard four-point backward 
difference expression for a second derivative. In defining the initial conditions z i n and 
Zi q, only the displacement quantities which are the first three elements of z need to be 
prescribed since the components in the last row in matrix D are all zero. 

Thus, by utilizing equations (10) to (19), the governing equations become 
Bl z l + A lZ2 = gi 

AlZi+l,* + Bizi t i + Cizi-i^ = gi (i = 2, 3, . . ., N - 1)\ (26) 

C N Z N-1 + b N z N = SN 
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where 


A l = 


O ( j 0 . H (A . Ac] 

*"/ 2~ 


Bl = 




gj - - S2 0 f O 
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n ; j l + h l / L 

%A t + -J-) * t 


C N " 


£ 2 t 


Jl h l 


LVT 


Ai 


2 / + 2 


SN " l ~L ' ^L f L 


2Ei 

A i = ir +F i 


-4Ei 


Bj = — — - + 2 AGi - 2 Aa t D 


2Ej 

C i = ”a~ "* F i 


A 


> (27) 




Si = 2 Aei + 2 ADlfyz ifl _ t + 


%Z-2 6z M-3 


a l z i,0 



Thus equations (26) are the complete set of equations required to solve for the 
vector z^ ^ at all stations i out to and including time step l . 


Step Loads 

Step loads or suddenly applied loads represent a discontinuity in the load history at 
a point in time k and, in essence, impart an acceleration to the shell. At this point in 
time k, 

z i,k 
z i,k 
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but 


2 i,k * (29) 

The superscripts - and + indicate quantities obtained by taking the limits of the 
quantities at a time k from below and above, respectively. The acceleration z* k is 
obtained from equation (10) with loads at k + imposed along with the deflection and 
velocity at k. Thus, the problem becomes essentially another initial-value problem 
with initial conditions given at time k. With z^, and z^ vectors computed, 

the problem is continued similarly to the sequence of equations (20) to (25). 

COMPUTER PROGRAM 


Program Organization 

The program developed in this paper takes the basic program in reference 1 and 
modifies that program in two major ways. First, the program HGS of reference 1 is 
modified to include a time-step cycling procedure to account for the inertia terms. The 
second major change is that the Fourier coefficients for 11 variables Ng, N^g, 

Q^, M^, Mg, M^g, U^, Ug, W, and contained in appendix A are summed for 

the series truncated at some value n and for each point in time l . The value n at 
which the Fourier series are truncated is supplied by the user. 

To accomplish the first modification the following subroutines were altered or 
added to program HGS of reference 1. 

Subroutine Remark 


MAIN 


INPUT 


CALZ 


MAIN structures the program so that HGS will calculate coeffi- 
cients of the 11 output variables at all interior stations i and 
boundary stations 0 and L for each station in time l and 
each Fourier variable n. It has been modified to include the 
inertial terms and the time-stepping cycles. 

INPUT is a user -prepared subroutine which supplies the shell 
geometry. In addition, the time increment e is supplied by 
this subroutine for the time-dependent problem. 

The subroutine CALZ stores or computes z^g, z^g, and z^g 
and is the only new subroutine added to HGS. 
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ABCG 


The subroutine ABCG sets the coefficients a, y, "6, a, 
and /3 in equation (19) as well as elements of the A, B, C, 
and g matrices in equation (26). Note that B and g 
include time -dependent terms. 

OUTPUT The subroutine OUTPUT controls program printing. In addition, 

it stores the vectors and z i,Z-2 for use in 

time step l + l computations. 


The flow chart of the present program follows. 


START 


Read 

INPUT 


Set Fourier 
variable n 


Calculate 
matrix D 
(eq. (10)) 


Set point in 
time, l 


Call CALZ 


Call HJF 
at i = 1 


Do to statement 
label 2, n = 0, 
NFOUR 


r Call INPUT: Set shell] 
^ geometry and e 


Do to statement 
label 19, l - 0, 
TIME 


If l = 0 
or 

If load discon- 
tinuous at 
time step, l 


Define H,J,f, 
equation (12) 


Call EFG 
at i = 2 


Call FORCE 
at i = 2 


Call ABCG 
at i = 2 


Call INIT 


Set spatial 
station i 


Call EFG 


Call FORCE 


Define E,F,G, 
equation (10) 


Define vector e, 
equation (10) 


Define A,B,G,g, 
equations (26) 


Reads Namelist 
FIRST (that is, 
boundary con- 
ditions at 
station 0) 


Do to statement 
label 3, i = 3, 
N - 1 
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Solve ABCG equa- 
tions by Gaussian 
elimination pro- 
cedure 


Reads Namelist 
LAST (that is, 
boundary con- 
dition at i = L) 



The second modification, that of summing the 11 variables at each time step, is 
accomplished by a followup program called SUMUP. Here SUMUP is actually a sepa- 
rate program that uses the same control cards as the HGS program. Thus, there are 
two programs but only one deck. The program HGS stores the information for all 
required time stations l and Fourier integers n of the 11 output variables on tape. 

The program SUMUP rewinds and reads the tape. Then SUMUP performs the series 
summation of the truncated series by using the Fourier coefficient taped from HGS. 

A listing of HGS and SUMUP is included in appendix D. Note that all subroutines 
which are not discussed in this report contain COMMENT cards in the listing explaining 
their functions. In addition, an asterisk in the right-hand column denotes all the state- 
ments that differ from the program presented in reference 1. The glossary of FORTRAN 
names of the variables is given in reference 1. 

Input Data 

The input data are those of reference 1 with some deletions and additions. These 
input data are now read into program HGS through the Namelists INPUTD, FIRST, and 
LAST. The variable type R stands for a real (floating point) value and I denotes an 
integer quantity. The complete list of Namelist INPUTD quantities and their definitions 
follow: 
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FORTRAN 


Name 

variable 

type 

Description 

NU 

R 

Poisson's ratio 

TKN 

R 

reference thickness, h Q 

CHAR 

R 

reference shell dimension, a 

EALSIG 

R 

Ea/cr 0 , thermal coefficient used only if thermal 
stresses are calculated 

IND5 

I 

an integer value of 0, 1, 2, or 3 

= 0, no poles; 

= 1, pole at £ = 0; 

= 2, poles at £ = 0 and | L 
= 3, pole at £ = | L 

NMAX 

I 

total number of meridional stations 

FREQ 

I 

integer which controls the frequency for printing 
numerical results. Results are printed at every 
FREQth station. 

NTIME 

I 

number of time steps of size e which will be taken 

JUMP 

I 

integer array (five elements) defining time stations 
at which a load is suddenly applied. If not 
required, then set JUMP(i) > NTIME. 

RUNTYPE 

I 

integer defining type of problem to be run: 

= 1, static case 

= 2, dynamic-response problem with either z^ q 
and/or z^ q specified by functions ZINIT 
and Z DINIT 

= 3, dynamic response problem 

NFOUR 

I 

last Fourier value to be summed 

KTHTIME 

I 

time-station interval at which coefficients of the 
11 variables from N = 0 to NFOUR will be 
printed out. If none other than at the zeroeth time 
station are desired, set KTHTIME > NTIME. 
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The namelists FIRST and LAST describe the boundary conditions at station 0 and 
L, respectively, and are required when IND5 of Namelist INPUTD fails to define a pole 
point at that station. These namelists define the elements of the matrices in equation (14). 
The input quantities of Namelist FIRST are 


Name 

Type 

Description 

OMEG1 

R 

a 4 x 4 array defining the 


matrix of equation (14) 

CAPL1 

R 

a 4 x 4 array defining the 

A 

matrix of equation (14) 

ELI 

R 

a 4 x 1 array defining the 

l 

vector of equation (14) 

Similarly, the Namelist LAST quantities are 



OMEGL 

R 

a 4 x 4 array defining the 

ft 

matrix of equation (14) 

CAPLL 

R 

a 4 x 4 array defining the 

A 

matrix of equation (14) 

ELL 

R 

a 4 x 1 array defining the 

l 

vector of equation (14) 


In addition to these input values there are various user -prepared subprograms. 
These subprograms, for the most part, are adequately described in reference 1. They 
include: 

Subroutine INPUT(NMAX) defining the arrays of p, y, u)g, c and 

the increments A and e 

FUNCTION HHT(K, DEL) defining h/h Q 

FUNCTION DHHT(K, DEL) defining (h/h 0 )' 

FUNCTION HRA(K, DEL) defining h/t 

FUNCTION DHRA(K, DEL) defining (h/t)' 

FUNCTION P(K, DEL) defining p( n ) 

FUNCTION PX(K, DEL) defining p^ n ) 

FUNCTION PT(K, DEL) defining p^ 

FUNCTION ZINIT(KK) defines z i Q 
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FUNCTION Z DINIT (KK) 


defines z 


i,0 


where 



KK = 1 corresponds to u^ 
KK = 2 corresponds to u^ 
KK = 3 corresponds to w 

and 

and 

and 

FUNCTION TEMP(K, DEL) 

defining T^ 


FUNCTION DELT(K, DEL) 

defining AT^ 


FUNCTION DTEMP(K, DEL) 

defining T^ 


FUNCTION DDELT(K, DEL) 

defining AT^ where h and t 


and face sheet thickness, respectively. 


The last four functions define the quantities and derivatives of the temperature 
equation 

T (n) _ t^(^) + AT^U) £ (30) 

where T^ n \|) is the Fourier coefficient of the temperature at the reference surface and 

ATj n \|) is the temperature difference between the inner and outer surfaces per unit 
thickness. 

The two additional functions, ZINIT and Z DINIT, have been added for use with the 
initial conditions required for RUNTYPE = 2. These functions define the initial 
replacements and velocities required by the initial conditions. 


Program Output 

The output comes from both programs HGS and SUMUP. The output from the pro- 
gram HGS consists of a complete list of the input data and of the shell geometry in tabular 
form (that is, column format) at all stations. This output is followed by tabular listing of 
the Fourier coefficients of the 11 variables in program HGS at the FREQth stations and 
at the KTHTIME cycle for each value of n. The output of program SUMUP then follows 
with the complete summation of the Fourier series to n equals NFOUR for the 11 non- 
dimensional variables at each KTHTIME cycle. The printout of SUMUP is at a value of 
8 where 


8 = AK v 

and AK is an INPUTD quantity. 


(31) 
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Estimation of Increment Size 


A discussion of the meridional increment size is contained in reference 1. The 
size of the time increment e from equations (19) to (24) also affects the results. Basi- 
cally, as the increment e tends toward infinity, the dynamic response is damped out. 
Thus, the increment e must be made sufficiently small. As shown in reference 4, the 
Houbolt method of initiating the problem and using backward differences is a stable 
method of numerical integration for transient problems. The selection can be made by 
comparing the results of increment size where increment is decreased by one-half until 
agreement between the results is obtained within suitable limits. 

EXAMPLE PROBLEMS 
Cylindrical Shell 

The purpose of this example is simply to demonstrate the accuracy of the program. 
The problem is that of a simply supported cylindrical shell loaded laterally with an axi- 
symmetric sinusoidal pressure oscillating with respect to time. The simply supported 
boundaries are free to displace in the u^ -direction at each end. The pressure q(£,t) 
becomes 

q(£,t) = p^ sin ^-cos 2nfl 
^L 

For this particular example, the frequency f and pressure amplitude p^ are taken 
to be 


f = 10 000 Hz 

p(0) = 1 

The geometric parameters of the shell are 
s/r = 10 
r/h = 100 
h/t = 2 
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Since at T = 0, there is an applied load on the structure, the initial displacements 
are assumed to be 


TT - 

U£ = u cos ~ 


*1 


U 0 = O 


W = w sin — 

where the amplitudes u and w are calculated from exact linear theory (ref. 3) to be 
u = 0.107388 

w = 1.11396 

Further, the initial velocity becomes 
& i,0 = 0 


for all i. 

In figure 3 the dynamic response of the deflection W at £, = 0.5 is compared 
with the exact solution by axisymmetric linear theory. (See ref. 3.) The curve in fig- 
ure 3 shows that this numerical solution and the exact solution are in close agreement. 
The numerical input for this problem was 

NMAX = 21 
v = 0.3 
AT = 0.1 sec 


p = 0.000259 

in^ 
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Figure Itynamic response of the deflection W at l L /2 
for cylindrical shell example. ' 


Conical Shell 

Description . - In order to demonstrate the data preparation required for using this 
program, a practical aeroshell problem is analyzed. The simply supported 120° trun- 
cated conical shell of sandwich construction studied in reference 8 has been selected. 

For this shell the loading is the aerodynamic pressure q acting laterally on the shell as 
it passes through the atmosphere. In addition, the shell axis has a small wobble or angle 
of attack \[/ which oscillates as a function of time, and thereby causes the loading to be 
applied asymmetrically. The purpose of such an analysis could be to ascertain whether 
the oscillations of the shell axis cause any stress buildup in the stress resultant Ng. It 
would be expected that any increase in Ng would have an important effect on shell 
instability. The shell cross section along with its physical dimensions is shown in fig- 
ure 4. The face sheets are made of aluminum and have a density p„ of 0.1 lb/in^ 

(2.8 Mg/m 3 ) and the core honeycomb has a density equal to 0.03p a . Thus, the average 
density of the shell p is given by 
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Figure 4.- Physical dimensions of the conical shell. 


P = P 


aVh 


— + o.03- — 2t 


(32) 


The dynamic forces are the Newtonian impact forces on the shell as it passes 
through the atmosphere with a small wobble \p and are defined to be 

p = -2q cos 2 /3 (33) 


where /3 is the angle between the normal to the shell reference surface and direction of 
the fluid flow and q is the aerodynamic pressure. Since <p is the colatitude angle 
(see fig. 1) of a point on the shell meridian, equation (33) becomes 


P = 



sin 2 ip sin 2 0 + cos 2 ip 


cos 2 <p 


2 


cos \p cos cp sin \p sin (p cos 6 



sin 2 i p 


sin^ cp 


cos* 4 6 


(34) 
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Thus from equation (34) and appendix B, the Fourier pressure coefficients become 
= -q^sin 2 0 sin 2 <fi + 2 cos 2 if/ cos 2 0^ 

p(!) = -4q^cos if/ cos 0 sin if/ sin 0^ \ 

p( 2 ) = -q^sin 2 0 sin 2 0 ) 

J 

Here the angle 0( t) is given by 


( 35 ) 


0 = 0 sin 277ft 


(36) 


where 0 is the amplitude of the oscillation and, for this example, is 5° or it/ 36. The 
frequency of oscillation f is 10 Hz. 

From equation (12) the simply supported boundary at 4=0 yields the conditions: 


0 

0 

0 

0 


0 

0 

0 

0 


0 

0 

0 

0 


r 

t 






♦ 


Jb L. 


0 

1 

0 

0 


0 

0 

1 

0 


0 

0 

0 

1 


r ^ 


u 0 

w 


r -\ 

0 


= 


m 
v v 


0 

v J 


(37) 


The edge at 4 = 4 l is pinned and restrained from horizontal displacement but 
allowed to displace in the vertical (axial) direction; thus, 
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(38) 


The initial conditions are taken as 


z i,0 = z i,0 “ 0 


(39) 
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Na melist INPUT . - In keeping with the definitions in section on input data, the 
Namelist INPUTD values become 


NU = 0.32 
TKN = 0.54 
CHAR = 90. 

EALSIG = 0 
IND5 = 0 
NFOUR = 2 
NMAX = 76 
FREQ = 2 

JUMP(l) = 82,83,84,85,86 
NTIME = 81 
RUNTYPE = 3 
KTHTIME = 82 
AK = 0. 

The boundary conditions required for Namelists FIRST and LAST have been given 
in equations (37) and (38) and correspond with equation (12). In Namelist format the 
values for Namelists FIRST and LAST are: 

OMEGl(l,l) = 16*0. 

CAPL1(1,1) = l.,4*0.,l.,4*0.,l.,4*0.,l. 

EL1(1) = 4*0. 

OMEGL(l,l) = 2*0. ,-0.5, 7*0. ,0.866, 5*0. 

CAPLL(1,1) = .866, 4*0., 1.0, 2*0. ,.5, 6*0. ,1.0 
ELL(l) = 4*0. 
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User -prepared subprog rams.- The pertinent statements for the conical geometry 
to be described in subroutine INPUT (see appendix D) are 

RI = 40.1 

RO = 90.0 

PI = 3.1415926535979 
ANG = PI/6 

DEL = (RO-RI) / (RO*COS(ANG) *FLOAT(NMAX-2)) 

R(NMAX) = 1.0 
R(I) = RI/RO 

DELR = (RO-RI) /FLOAT(NMAX-2) 

R(2) = R(l) + DELR/2 
NM1 = NMAX -1 
DOlI = 3, NMI 

1 R(I) = R(I - 1) + DELR 
D02I = 1,NMAX 

GAM (I) = COS(ANG)/R(I) 

OMT(I) = SIN(ANG)/R(I) 

OMXI(I) = 0. 

2 DEOMX(I) = 0. 

C IF RUNTYPE = 2 OR 3 THEN SET EE AND RHO TO FIND TIME INCREMENT, EPS. 
C HERE EE IS THE MODULUS OF ELASTICITY AND RHO IS DENSITY (LB/IN**3) 
GACC = 386.088527 
EE = 10500000. 

RHOA = 0.1 
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RHO = RHOA*2.75/27. 

SS = EE /RHO 

EORHO = SS*GACC 

EPS = SQRT(EORHO/CHAR**2) *TIM 

A complete listing of the program is contained in appendix D. The statements 
required for the function subprograms become 

HHT = 1.0 

DHHT = 0. 

HRA = 27. 

DHRA = 0. 

PX = 0. 

PT = 0. 


TEMP = 0. 

DELT = 0. 

DTEMP = 0. 

DDELT = 0. 

ZINIT = 0. 

Z DINIT = 0. 

It should be noted that the left-hand side of these statements defines the function subpro- 
gram in which the statement is used. The statements required for function P(K) are 
more involved. By setting q in equation (35) equal to unity and recalling equation (36), 
the required statements become 

PI = 3.14159265358979 

CPS = 10. 

TIM = 0.05/CPS 
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TAU = FLOAT (JTIME)*TIM 


AO = 5. 


A = AO*PI/180. 

FEE = PI/6. 

PSI = A *SIN (CPS *T AU *2 . *PI) 

Q = l./LAM 
SA = SIN(PSI) 

CA = COS(PSI) 

SF = SIN(FEE) 

CF = COS(FEE) 

PO = -Q*(SA*SA*SF*SF + 2. *CA*CA*CF*CF) 

PI = -4.*Q*CA*SA*CF*SF 
P2= -Q*SA*SA*SF*SF 
IF(N.EQ.0.)P = PO 
IF(N.EQ.1.)P = PI 
IF(N.EQ.2.)P = P2 

Output .- The results generated by the program HGS are summed by the followup 
program SUMUP. The tabular results at 0 = 0 (that is, AK = 0 in Namelist INPUTD) 
are printed for the 11 variables and are shown for time t = 0. It should be noted that the 
output is nondimensional. In other words, the multiplicative constants involving a 0 , h 0 , 
a, and E 0 in front of the summation signs in appendix B are not included for the output 
results. In order to get dimensional results, these nondimensional results must be multi- 
plied by the appropriate constants of appendix B. 
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For this particular example, the maximum stress resultant value of N$ which 
occurs at | « 0.70 £l is plotted in figure 5. Also included are the results when the 
oscillation frequency f is 1 Hz or 100 Hz. In figure 5 the time axis is made nondimen- 
sional by T, the time required for one cycle of oscillation (T = 1/f), for ease in com- 
paring different oscillation frequencies. The time increment e taken for this example 
was 0.05T. The difference in the results for frequencies of 1 Hz and 10 Hz is negligible 
whereas for the 100 -Hz case, there is substantial increase in compressive stress. The 
time-response values at 1 Hz and 10 Hz settle down almost immediately to a harmonic 
steady-state response pattern whereas the 100-Hz curve is still nonharmonic after four 
cycles because of the larger change in z^ q. The change in pressure is occurring so 
rapidly in this case that the rate of change in loading is approximating a step load. 



Figure 5*- Maximum stress resultant (£ ** O.JO and 0 = 0°) as a 
function of time for various oscillation frequencies. 


CONCLUDING REMARKS 

This report describes the development of a computer program for the linear asym- 
metric bending behavior of a statically or dynamically loaded elastic thin shell of revolu- 
tion subjected to either thermal or mechanical loads. The program is an extension of the 
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static analysis described in NASA TN D-3926 to include dynamic loads and summation of 
Fourier components. These changes required the addition of the dynamic terms in the 
equations as well as the restructuring of the program to include the initial conditions and 
backward numerical integration of time derivatives. Two examples, demonstrating the 
flexibility of the program as well as the data preparation required, have also been 
included. The first example, response of a thin cylindrical shell to an oscillating pres- 
sure, demonstrates the accuracy of the program for dynamic -response analysis. The 
second example is used to demonstrate the preparation of user-supplied data for a prac- 
tical analysis. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., September 22, 1970. 
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APPENDIX A 


CONVERSION OF U.S. CUSTOMARY UNITS TO SI UNITS 

The units used for the physical quantities in this paper are given both in the U.S. 
Customary Units and in the International System of Units (SI). Factors relating the two 
systems are given in reference 7 and those used herein are given in the following table: 


Physical quantity 

U.S. Customary 
Unit 

Conversion factor 
(a) 

SI Unit 
(b) 

Length 

in. 

0.0254 

meters (m) 

Density 

lbm /in3 

27.68 X 103 

kilograms/meter 3 (kg/m3) 

Modulus, elastic 

psi = lbf/in^ 

6895 

newtons/meter2 (N/m^) 


a Multiply value given in U.S. Customary Unit by conversion factor to obtain 


equivalent value in SI Unit. 

bPrefixes to indicate multiple of units are as follows: 


Prefix 

Multiple 

micro ( jll ) 

CO 

J 

o 

r-H 

milli (m) 

10-3 

kilo (k) 

103 

giga (G) 

10 9 
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APPENDIX B 


FOURIER SERIES EXPANSIONS 

The Fourier series expansion of the dependent variables in the circumferential 
direction is presented in this appendix. The constant terms to the left of the summation 
symbol are required to nondimensionalize the series coefficients in a consistent manner. 
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APPENDIX C 


DEFINITION OF MATRICES 

As shown in reference 2, the nonzero elements of the matrices E, 

are: 


E 11 *b 

- b ( 1 ~ v ) + >; 2 d(l - ^)(3a)g - 


e 22 

e 23 


A 2 d (l - ^)(3o)g - co^n 

2p 


E 32 - e 23 


e 33 = A 2 d(l - v) 


+ (1 + i/)y2 

P l 


E34 = A 2 


e 43 = ~d 


Fji = yb + b' 




e 13 = b(w| + vuiQj + A 2 d(I - v) 
e 14 = 


/ 2 \ 

(1 + ^y 2 ^ + J(3co^ - c o e ) 


G, and e 
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APPENDIX C - Continued 

F 21 = _F 12 

2 

f 22 = + b ’) “ X 3 ~ ~ ( 3aJ e " - rfaz - 3w e )J 

X 2 d’(l - v)[ \2 

+ 8 ( 3 "« ■ "«) 

x 2 d(i - y) n r . . , 0 , <r\ * 2 d’(i - y )( 3w 0 - w ^) n 

F 23 = 2 p L 2(1 + l,)r ° Je ~ + M"* ~ Wfl )J + 2p 
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f 32 - - *»«<« - 2 ") - “{'] + V) > e - » { ) 

2 “ g 

F 33 = - v) (1 + v)(2y(jL> ^OJ q + y^j + ^ + X 2 d’(l - z-0 (1 + f)y 2 + 

f 34 = A 2 y(2 - v) 

F44 = da>^ 


F 4 o = -d^y 


, 9 (1 - f)bn 2 

Gn = vb y - PbcogW^ - by^ - 

9 9 2 ( 3w £ ~ ^q) 2 ^ 

- X 2 d(l - v) (1 + v)y 2 w? + ^—2 

5 8p 2 
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12 p 2 p 


X 2 d(l - v)yn 
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( 3w 0 - w $) 
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+ (1 + V)(x>£COq 


n u T 1 ! \] , u if \ X 2 d(l - z^yn 2 3co £ x 
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G 14 = * 2 U " v)yu)£ 
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APPENDIX C - Continued 
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J 32 


= bn(o>0 + vco£ X 2 d (i _ v)n 


2p 


n _ 

2(1 + ^)(o)|O)0 - y 2 o>£ + 2y 


‘co 


n 2 co 0 \ 


P 2 


yo)| + 3y 2 (a>0 - co ^ + cog(Og(3co 0 - 00 A 

~ d ' 2p~ ' "^ )n j^ (1 + v)yu> e + y ( 3a, 0 - 


^33 + 2pC0^U)q + O)0 2 ^ + - — - — - — (1 + v ) (co^cog - ~ + 2y 2 ^ 


2 (y 2 + CO j: CO qJ 


A 2 d'(l - j))n 2 


(3 + v)y 



where 
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G 34 = 


->2 


(1 - 


G41 

= d|o)| + 

pyoj^j 



dvnw Q 



g 42 

_ (7 

P 



g 43 

dm 2 

" P 2 



g 44 

= -1 



e l = 

- P£ + tj 1 

- X 2 (l - 

U)yC0^Mrjr 

e 2 = 

-n - 

r - a 2 (i 

- v)^<x> e M T 

e 3 = 

-p - (■»{ 

+ 0)^t T 

- X 2 (l - v)yMr^ 

e 4 = 

M'p 








[M'p 


E 0 h 0 (l - v 2 ) 

y E£ 2 d? 


~ E 0 h 0 3 (l - v 2 ) 

j* EaT^ d? 
tT ^oV 1 ‘ v ) 
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Mrp = 



^Eq-tW d£ 
a 0 h 0 3(i - v) 


The nonzero components of the H,J,F matrix (ref. 2) are 


H n = b 


h 23 - X2d(1 2p : t ' )n (3“e - "|) 
H32 = ( 3 “e - "i) 


H 33 = X2d(l - v) 


^ + <1 + vrf 

p* 


H34 = A 2 

H 43 = -1 

J-q = vyb 

T ^nb 
J 12 = — 


J13 = b ( w £ + vw e) 


*21 


b(l - v)n dA 2 (l - v)n 


2 p 


8 p 


( 2a, £ - co 0 )(3u> 0 - co^ 
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J 22 = -r 11 22 


J 23 = 23 


J3I = -A 2 d(l - v) (1 + ^)y 2 co£ + - <*>e) 


J 32 = ” 


X 2 d(l - ^)yn 


— [3a>0 - W| + 2(1 + 


J 33 = -A 2 d(l - v )(3 + v)^L 


J34 = X 2 (l - v)y 


J 41 = 

Fl = -t T 


F 3 — X 2 y(l - I'') XI rp 
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APPENDIX D 


PROGRAM LISTING 

The complete listing of the program with both the HGS and SUMUP programs follow. 
The asterisks on the right-hand edge indicate those statements added to or modified from 
those given in the HGS program of reference 1. 
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PROGRAM HGS C INPUT tOUTPUT ,TAPE5=INPUT * T APE6=0UTPUT , T AP £ 1 A ) 

C MAIN PROGRAM 

C THIS PROGRAM CONSISTS OF THE MAIN PROGRAM TOGETHER WITH THE FOLLOWING 

C SUBROUTINES 

C MAT I NV 

C EFG 

C OUTPUT 

C STRESS 

C EQ73 

C INI T 

C FINAL 

C FORCE 

C PANDX 

C ABCG 

C HFJ 

C KLT 

C BOB 

C POLE 

C IN ADDITION TO THE ABOVE THE USER MUST SUPPLY THE FOLLOWING SLB 

C ROUTINES ANE FUNCTIONS. 

C SUBROUTINE INPLT — CALCULATES THE SHELL GEOMETRY. 

C FUNCTION P-SPECIFIES THE NORMAL PRESSUSIE DISTRIBUTION. 

C FUNCTION PX-SPECIFIES THE MERIDIONAL PRESSURE DISTRIBUTION. 

C FUNCTION PI -S FECI FI ES THE CIRCUMFERENTIAL PRESSURE DISTRIBUTION 

C FUNCTION TE^P-SPFCIFIES ThE MERIDIONAL TEMPERATURE DISTRIBUTION 

C FUNCTION GTEMP- SPECIFIES THE DERIVATIVE OF THE TEMPERATURE 

C FUNCTION CELT- SPECIFIES THE DISTRIBUTION OF THE TEMPERATURE 

C VARIATION THROUGH THE THICKNESS. 

C FUNCTION CCELT- SPECIFIES RATE OF CHANGE OF DELT 

C FUNCTION HHT-SPECIF IES THE SHELL THICKNESS DISTRIBUTION 

C FUNCTION CFFT—SPECIFI ES THE DERIVATIVE OF THE THICKNESS. 

C FUNCTION HPA-SPEC IF IES THE DISTRIBUTION UF THE RATIO OF THE 

C TOTAL THICKNESS - TO -COVER PLATE THICKNESS. 

C FUNCTIGN DHP A- THE CERIVAflVE OF THE ABOVE RATIO. 


INTEGER FREGfCYNRtSPfBCUNDRY , RUNTYPt * 

REAL NU,LAM,N,JAY 

COMMON 

3/BL3/NU,LAM,N,EALS IGt CHAR, DEL 
6/OLk/A(4»A ),e(A,A),C(4,A) f S MAGI 4) 

7/BL 7/PEC { A f A , 1C2) ,X< 4, 1C2 ) * 

9/BL9/D (A, A ) ,ZC0T ( A, 102 ) , ZDDO T ( A , 10 2 ) , EPS * 

A/BLiO/JTIME * 

COMMON /BL1 1/ZSAVfc (A, A, 102) * 

C/BL12/INITZ * 

o/BL 14/AC, 7 I M , CP S * 

G/BLib/TAU * 

H /BLI7/KTHT IME * 

K/BL18/RLNTYPE , JUMP(5) * 

DIMENSION 2(4,102) ,IPIVCT(A) ,INDEX(4,2I * 

EQUIVALENCE (X ( 1, 1 ) ,Z( 1, IH * 

NAMFLlST/IAFLTC/NLfTKN.CHAR, EALS IG, INC 3, * 

1NFUUR, NM A X , FR EC, JUMP, NT IME , RUNT YPE , KTHT IME , AK * 

C RUNTYPE ^ I STATIC CASE * 

C RUNTYPE = 2 DYNAMIC RESPONSE — INITIAL DEFLECTIONS * 

C RUN1YPE = 3 DYNAMIC RESPONSE' — INITIAL LOADS * 

12 KE AC ( 5,INFLTC > * 

PRINT I N PUT C * 

RT=RUNTYPE * 

READ( 5,31 ) 

31 FORMAT (72H 

1 ) 

WR I TE< 6,31 ) 

JT I ME=0 * 

N = 0 • 

L AM- TKN/CHAP 
K - 1 
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2U FORMAT ( 1 rl l * 1 *3 h TIMF I TF RAT I ONI 5/1 OH REAL T I MF 3X 1 1 5. 8 , 5 X*P ( 1 ) = * * 

1 h 1 3 - 7 I * 

3u CONTINUE * 

IF ( jr IPt • NE . C ) GO TO 17 * 

IF (KUMYPE.EC*! ) GC TG 16 * 

IF (kUMYPE.EC.3) GU TO 16 * 

CALL CaLZ (WAX, INL5) * 

NPkJNT=2 

CALL OUTPUT ( FR EQ , NMAX , DEL , N PR I NT ) 

17 CONTINUE 

IF ( JTIPb.NE. JLPP< JJJ) ) GU TO 16 * 

CALL CALZ (WAX,INC5) * 

JJJ=JJJ+1 * 

Lo NPRINT=3 * 

IF (RUMYPE . tC.2 . ANC. JTIMP.tQ.O) GO TO 19 * 


CALL HFJ ( 1 , JNC5,WAX,2.) 

CALL tFG(2,INCb,NPAX) 

CALL FCKCE ( 2 , IND5 ,WAX ) 

K = Z 

CALL ABCG ( K ) 

CALL IMT ( INCb,HCUNDKY ) 
l\MAXi=NPAX-l 
U 0 6 K = 3 , N P A X 1 
CALL EPC(K , INC5, WAX) 

CALL FCRCE(K, INC6,WAX) 

CALL ABCG ( K ) 
o CALL PANCX ( K ) 

CALL HF J( NPAX , IN05,KMAX f 2. ) 

CALL FINAL (WAX, I N C b , P CUN C R Y ) 

UL 7 L=2»NPAX1 
K=NPAX + 1-L 
/ CALL tUi(k) 

IM INU*>.eC.C.CR.INC5.EC.3) GU TO 14 
CALL L C 7 3 ( 1 ) 

CL TC lb 

14 CALL EFL( 2 , INC5,NMAX) 

LALL FORCE (2, IN Lb, WAX) 

K=2 

CALL A B C C (K ) 

Uu d I = L , 4 
b 1 = G . 
b2 = 0. 

UN 4 J= 1 , 4 

Sl= Sl + M I , J)*Z ( J,3 ) 

4 S 2= ;>2 + p ( I ,J)*Z(J,2) 
b SKAG1 I) = SPAC ( I J-S1-S2 

CALL MA1I NV(C ,4,SPAC, 1 ,UCTFRM, I PI VCT, INCtX,4, I SCALE) 
UC 10 1=1,4 

1 u L (I ,1) = $P A C ( 1 ) 
lb CALL STRESS (FPEC, WAX, IN Ob) 

CALL OUTPUT (FREO, WAX,L'FL, NPRINT > 


IF (KUNTYPE • N L • 3 ) GO TO 21 * 

IF (jriNt . N F . JUPP (JJJ)) GO TO 21 * 

INI TZ= 1 * 

LALL CALZ(NPAX, INCb) * 

NPKlNT=2 * 

CALL OUTPUT (FREC, WAX, OE L » NPRINT) * 

JJJ=JJJ+1 * 

21 CONTINUE: * 

bOUNDKY=l * 

IF ((KUNTYFf .EC. 1) . A MJ . (N . EQ . FLU AT ( NFOUR ) ) ) STOP * 

IF (RUMYPE.MT. 1) RUN T YP E = 3 * 

IV CLNTINLE * 

IF ( N .EG. FLCATI NFCUK) ) STOP * 

2 CONTINUE * 

CL» TG 12 
t. N 0 
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FHH J = H FT ( KfDtL ) 

FCHHT= DHFT ( K * CFL ) 

FHkA= FRA(K,CEL) 

FDHkA=CFRA!K, CEL I 
w=P!l) 

4 FORMAT ( 5X * *F/FU= * F 1 2 . 5 f 5 X * E ( H /HO ) /DS=*E 1 2. 5 f 5X*H/T = *E 12 • 5 , 5X 
1*D(H/T)/CS=*F12*5) 

WR I T t ( 6 f 4 ) F F F 1 , FDHHT t F H R A , FCHRA 

5 FORMAT !5X*TFIS RUM WILL SUV FOURIER TERMS FROM N=0 TO *14) 

WRITE! 6 ,5 ) NFCLR 

WRITE! 6,81 ) TIM, CPS 

81 FORMAT (5X*1 I ME I NCR EM t NT = * 1 1 4 . 7 , b X *CP S = *E 14 • 7 ) 

REWIND 14 

CALI RECOOT ( 14, 1, C, FREQ, N M A X , NTIME, NFOUR , AK, TIM, RUNTYPE) 
NF 1 = NFCUR+1 
DO 2 IN =1, N F 1 
IM=I N-l 
R UN T YP E = R 7 

IF (RUN TYPE. CC*1 ) WRIT E(6, 3 2) 

32 FORMAT(/16F STATIC PRUBLtM) 

IF (RUN TYPE. EC. 2) WPITE(6,33) 

33 FORMAT ( / 5 6 F DYNAMIC RESPONSE PROBLEM WITH INITIAL DEFLECT I C N S GIVE 
IN) 

IE (KLMYPt.bC.3i WRITE (o, 34) 

J4 FORMAT !/75F DYNAMIC RESPONSE PROBLEM - UEF L ECT I CN S CALCULATED FROM 

1 INITIAL LCACS GIVEN) 

DU j 1=1,4 

DCJ 3 J = 1 , 4 
D ( I , J ) = C . 

DC ltl) = l. 

D<2,2)=1. 

D ( 3 , 3 ) = 1 . 

DO 1U1 K = 1 , NM A X 
DC 101 J=1 , 4 

Ptt(J,3*K)=C. 

ZSAVE! JtltK)=C. 

ZSAVF ( J ,2 , K ) = C . 

ZSAVE ( J, J, K ) = C « 

ZS % VE t J ,4 , K ) = C . 

Z DU T ( J , K ) = C . 

101 ZDDOTt J,K )=C. 

CALL INPUT (NMAX, NTIME) 
wRIIE ( 6 , 6 C ) 

60 FORMAT t/llF INPLT DATA) 

hK i TE( 6,164 ) NL»TKN',CHAP»N, t ALS I G , NMA X , FR EQ 
164 FORMAT ( /I0>26FPCI SSC>iS RATIO ( NU ) =£16.8/ 

15X ilHREFtRFNCE THICKNESS (CKiM) =£16.8/ 

2 8X2 8HR EFERENCE LENGTH (CHAR) =£16.8/ 

31 1X25 H FOUR IbR INDEX (N) =£16.8/ 

41X36HTEMPEPATLPE COEFFICIENT (FALSIo) =£16.8/ 

56X3 OHN L m 8 t P OF STATIJNS ( NMAX ) =116/ 

6 6X3 Oh PRINT ING FREGLENCY (FREQ) =116) 

IF (kUNTYPE.NE.l) WRITE(6,35) NTIME 
36 FORMAT (4X 32FNLMBER OF TIME STEPS (NTIME) =116) 

BUUNOR Y = 0 

IF ( IN .NE. U BCUNURY =1 
N P R I N T = 1 

IF (N .EG. O.KALL CUTPUT(FrEQ, NMAX, DEL, NPPINT) 

J J J = 1 

DU IH I 1= L , N T IMF 
J T I ME = I I— 1 
I N I T 7 = 0 
0 = P ( 1 ) 

IF (MUC ( J T IME ,KTFT I ME ) .NE . 0) GO TU 80 
WK I TF ( 6 , 2 G ) JlIMF ,TAU, Q 
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S L *3 k u U T IKE Lfill (NPAX^IND5) * 

f\ h AL \L»i\»LAP * 

CCMVufM * 

1/iiL 1/K( IG2 ) *GAM 1C2) tUMT<lC2 ) ,QMXI ( 102 > , DhCMXf 102 ) * 

3/t3L3/NU,L AP f N, EALSZG,CEIAR, DEL 
6/HLO/A ( A ,4 ) f B (4,4 ) ,C( 4 f 4) v SJ*AG(4) 

7/l3L7/PEE(4,4 # lC2> »X(4,IC2) * 

<J/riL9/L)(4t4 ) »ZCCT t 4 ,10 2 ) t ZDCOT { 4 , 102 ) , EPS * 

A/3LL0/JIIPE * 

3/BL 1 1/ /SAVE ( 4,4,10?) * 

C/HL12/IMTZ * 

.< / *3 1 1 3/HUN TYPE . JUPP(5) * 

IKTLGEP KLMYFE * 

IE ( RUN TYPE JF. ?) GJ TJ U * 

C Ptt ( 1 , i ,K) =U X I * 

C P fc E ( 2 i 3 , K ) = L THETA * 

C Pth ( J » 3 » K )" w * 

C ZOuT ( l,K)=b XI r.CT * 

C ZUuTt 2,K)=b TEt IA HOT * 

C ZOuT ( 3 * K ) ■= V, CGT * 

UU I K=1,KPAX * 

DU 1 J = 1 , 3 * 

PfLlJt 3*K) = HMT(J,K) * 

1 ZiKlK J f K) = ZDlMT( J,K) * 

KM A X 1 - N y A X - l * 

K<U 1 = ( Pt.E ( 1 , 3 , 1) + PEE( l ,3 ,2 M /2. * 

PR )2= (PEE (2*3 f 1 ) 4 P F F ( ’» 3»Z ) ) / 2 • * 

Prt'M=(PfE(3,3»L)+PEE( 3,3,2) ) /2 • * 

PPL t* ( PEE ( 1,3 » N M A X ) +P E E ( 1 , i ,KMAX1 ) ) /2. * 

Pr<L2=(PtH (2,3,ENAX)4pEE(2,3,NMAXl) )/2. * 

PKl. 3= ( PL t i IS 3 ,KMAX ) 4PHF ( i , 3 * KMAX1H /2 . * 

C'W-CUL A1 I u'\ Lf P X I * 

DL«:K = lfK*AX * 

IF IK.FUII GL TG 3 * 

IE (K.EC.NEAX) GO TO 4 * 

wL)P = ( P tH 3 ,K4 1 )-?.*PEE ( 3, 3 ,K ) 4PEE ( 3, 3, K-l ))/DEL**2 * 

aP= (Hf E ( 3 f ^ , K + 1 ) — P F F ( 3, 3» K— 1 ) ) / ( 2 * 3(5 Dc L ) * 

JAi P = ( PN ( 1 f 3 f K + 1 )-PFf: ( If 3 , K-l ) )/( 2 .*uru * 

oC I u *3 * 

^ *P = (Pt t ( 3 ,3 ,2 )-PT F i 3, 3, 1 ) )/GEL * 

J X 1 r •= ( PF fc( l, a ,2) — Pff{ 1*3,1 ) )/ObL * 

*mp= < J .**H fr t (3, J, 1 > - 7* *P£ E ( if 3»2 J + 5 .*PEE ( 3, i» 3)-PCE( 3, 3,4)) / * 
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1(2.*DEL**2) * 

CALL BOB (K,DEL»NU,BLK1,BLK2,DM,BLK3) * 

PEE 14,3,1 ) = DH*( — WDP + QMXK 1 ) *UXI P+ ( DEOMX { 1 ) +NU*G AM ( 1 * ♦QMXI ( 11) * 

1*PRG1-NU*GAM( 1 )*WP+NU*(N/R( 1 ) ) **2*PRO 3 + NU* ( N/R ( 1 ) ) *QMT ( 1 ) *PR02 ) * 

PEE(4,3,1)=C.C * 

GO TO 2 * 

4 WP=(PEE(3,3,NMAX)-PEE(3,3,NMAX-1))/QEL * 

UXIP=(PEE( 1,3, NMAX )-PEE( 1,3, NMAX-1) ) /DEL * 

WDP=(3.*PEE(3,3,NXAX)-7.*PEE(3,3,NMAX-1)+5.*PEE (3,3,NMAX-2) * 

1 -PEE< 3,3,NMAX-3) )/(2.*DEL**2> * 

CALL BOB (K,OEL,NU,BLK1,BLK2,OM,BLK3) * 

PEE (4, 3, NMAX ) = DM* (-WDP+OMX I ( NMAX) *UXIP+ { DEOMX < NMAX ) +NU*GAM ( NMAX) * 

l*OMXI( NMAX) )*PRL1-NU*GAM(NMAX)*WP+NU*(N/R(NMAX) ) **2*PRL3+NU* * 

2 (N/R (NMAX) ) *OMT (NMAX) + PRL2 ) * 

PEE (4,3,NMAX)=C.O * 

GO TO 2 * 

^ CALL BCB (K,DEL,NU,ELK1,BLK2,DM,BLK3) * 

PEE( 4,3, K)=DMM-WDP+OMXI(K)*UXIP+( DEOMX ( K ) +NU*G AM ( K ) )*PEE( 1,3,K) * 

1— NU*GAM (K )*WP+NU* <N/R(K >)**2*PEE(3, 3,K ) +NU* ( N/R ( K ) ) *OM T ( K ) * 

2*PE E ( 2 , 3 , K ) ) * 

2 CONTINUE * 

PR04 = P E E(4,3,l) * 

PRL4=PEE(4*3,NMAX) * 

PEE(4,3,1 >=2.*FR04-PE£i4,3, 2) * 

PEE(4,3,NMAX)=2.*PRL4-PEE(4,3,NMAX1 ) * 

12 NMAX 1= NMAX- I * 

IF (RUNTYPE .EC. 2) GO TO 21 * 

DO 11 K=2,NMAX1 * 

DO 11 J=l,4 * 

11 ZDOT( J,K)=< ll.*ZSAVE( J , 4 , K ) - 18 . *ZS AVE ( J , 3 , K ) +9.*Z SA V E ( J , 2 , K >-2 . ♦ * 

1ZSAVE(J,1,K) )/(6.*EPS) * 

21 DO 9 K = 2 , A M AX 1 * 

CALL EFC (K,INC5,NMAX) * 

CALL FORCE ( K , I ND5 , NM AX ) * 

CALL A8CG ( K ) * 

DO 10 1=1,4 * 

SUM=0. 0 * 

DO 20 J= 1,4 * 

20 SUM=SUM+A(I,J)*PEE(J,3,K+1)+E(I,J)*PEE(J,3,K ) +C ( I , J ) *PEE ( J , 3 , K- 1 ) * 

10 ZDOOT ( I,K)=( SUM— SHAG ( I ) )/( 2.*0EL) * 

ZD DOT ( 4 , K ) = 0.C * 

9 CONTINUE * 

RETURN * 

END * 
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SUBROUTINE F F J < K , I ND5 , N MAX , Y AH ) 

C SUBROUTINE FFJ THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE H,F, 

C ANO J MATRICES, AS CEFINEO IN APPENDIX A OF REFER ENC E { 1 ) , A T THE STATION 
C SPEC I F I ED BY TFE INDEX K. 

COMMON 

1/BL1/RC102) »GAM( 1 C2 I , OMT ( 1 02 ) , OMXI ( 102 ) , DEOMX I 1 02 ) 

3/BL3/NU, LAM, N,EALSIG, CHAR, DEL 
5/BL3/H(4,4 ),FF (4) , JAY I 4,4) 

REAL NU,LAMfN, JAY, L2 

CALL BC6( K,DEL,NU,B,DB,D,DD) 

L 2= LAM**2 
Dl=(l.-NU) 

QX = C MX I (K ) 

REG=0 • 

I F ( YAH. EC. 2, ) REG = 1 • 

E AL=E A L S I G 
T = T EMP ( K , C EL ) 

DLT = DELT ( K , C E L ) 

H1=HHT(K,DEL> 

HRB=HR A ( K , L EL ) 

FF ( 1 )=-2.*F 1*EAL/ <D1*HRB)*T 
FF ( 2 )= 0 • 

FF(3)=L2*CFAR*EAL*CLT*Hl**3/3.*(L.5/HRB-3./HR6**2+2./HRB**3) 

F F ( 4 )- 0 . 

I F ( INU5 )I, 1,2 

2 IF ( ( ( I ND5-2 ) .LE .0 ) .AND. IK.EQ.1) ) GO TO 8 
IF ( { ( I ND5-2 ) .GT .0 ) .AND. (K.EQ.NMAX ) ) GO TO 8 
1 GA=GAM(K) 

FF ( 3 ) = F F ( 3 ) *G A 
R A- R ( K ) 

UT= OMT IK) 

ENR=N/R A 

UXT=3.*CMXI ( K )~CMT { K ) 

UTX=3 . *CMT ( K ) — CMX I ( K) 

DL=D*L2*D1*ENF 
H ( I , I ) - B 
H 1 1 , 2 ) = C . 

H( 1 ,3) =C. 

H ( 1 , 4 ) = C • 

H(2,1)=C. 

H(2 ,2 ) =B*D 1/2.+ L2*D*D 1/8. *0TX**2*REG 
H(2,3)=DL/2.*CTX*REG 
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H(2,4)=C. 

H ( 3 , 1 ) = C • 

H ( 3,2)=0L*CT**YAH/4. 

ENR2=ENR**2 

H( 3 ,3) =L2*C*Dl*(YAH*ENR2+( L • +NU ) *GA**2 ) 

G A2— GA**2 
H ( 3 , 4 ) = L 2 
H(4,1I=C. 

H(4,2)=C. 

H ( 4 ,3 J - — 1 • 

H(4,4)=C. 

JAY (l f 1 ) = M*GA*R 
J AY ( 1, 2) =NU*B*ENR 
J AY ( 1 « 3)=B*(CX+NU*0T) 

JAY ( 1 1 4 )= C • 

J AY ( 2 ? l)=-B*Cl*ENR/2.-DL/8„*CXT*0TX*REG 
JAY (2,2) =-G A*E (2,2) 

JAY( 2, 3 ) — -GA*H(2, 3) 

JAY(2,4)=C. 

JAY(3,1)=-L2*C*C1*(( L .+NU) *GA2*CX+ENR2/4.*OXT*YAH) 

J AY ( 3 , 2 )=— GA*DL/2**( 2-*OT*( 1 . +N U ) +UT X / 2 . *Y AH ) 

JAY (3,3 )=- L2*C*Ci*< 1 . + NU+Y AH ) *GA*ENR2 
JAY( 3, 4)=L2*D1*GA 
J AY ( 4 , 1 ) = C X 
J A Y ( 4 , 2 ) = C * 

JAY (4, 3 ) = C • 

JAY (4,4 |=C. 

GO TO 5 
t> DO 61= 1,4 
DO 6J=1,4 
H( 1 , J) = 0. 

6 JAY ( I , J) = 0. 

H( 1, l)=B*< l.+NU) 

H(1 ,2) = N * B * N U 
H(2,I)=H*D1*( — N ) / 2 • 

AW6 = D1* ( - 3 * * ( 1«+NU)+N**2*( 3 .+NU>) 

Ab = Ai*B*C 

L 1=AW/ (-2.+NU* (N**2-l . ) ) 

ATHETA=1. 5<C*MCX*D1*( 3.+NU ) 

A X I = 1 • 5$D*CX*C1*(2**( 1 • + NU ) — N ) 

H(3,4)=L2* (?.-NU+2.*AWB ) 

H(4,3)=-l. 

JAY(1,3)=B*< 1 • +MJ )*CX 
J A Y ( 4 , 1) = C> 

DH= CHH T ( N , C t L ) 

DTOH=DHRA ( ,CE L ) 

OU|_T=DDfcLT(R,CfcL) 

DTM CM = C H AR *E ALS I G/3 •/□!*( ( 1 . 5/HRR- 3 . / HR 3**2 + 2 . /HR B** 3 ) * ( H 1 **3*DDL T 
l+3.*OH*hl**2*CLT) + DLT*H l**3*DTOH/HRB**2*(-l . 5 + 6 . /HRB-6 . / 

2HRB **2 ) ) 

EEU)=DTNCM ( 2 . * A k B / ( -2 . + NU * ( N* *2- 1 . ) ) +2.-NU) 
b RETURN 
END 



APPENDIX D - Continued 


SUBROUTINE OUTPUT (FREC,NMAX,QEL,NPRINT) 

C SUBROUTINE OUTPUT THIS SUBROUTINE CONTROLS PROGRAM PRINTING AND 

C PUNCHING. GEOMETRIC DATA IS PRINTED IF I NO 1 IS NOT EQUAL TO ZERO. ANY 
C OR ALL OF TEE ELHVFN OUTPUT QUANTITIES CAN eE PUNCHEO BY SUITABLE 
C SPtCIF ICAT ION CF TFL FIELDS OF THE NOL CARD 
COMMON 

1/bL 1/R ( 1C2 ) *GAM( 1 C2) ,0MT<1C2) ,UMXI I 102) ,DECMX( 102 ) 

7/BL 7/ PEE {4 ,4, 1C2) , X(4 f 1C2) 

9/6L9/D (4,4) ,ZCCT(4,102 > , ZDDCT ( 4 , 10 2 ) , E PS 
A/BL10/JT [ME 
6/bLll/ZSAVE(4 i 4 * 1 C2 I 
H /BL17/KTM IME 

DIMENSION J J( 1 1 ) , KK ( ll)»ESS(102 ), YORDt 102) 

EQUIVALENCF ( X ( I , 1 ) f E SS ( 1 > ) , ( X (1 , 20» YURD ( 1 ) ) 

INTEGER FR EC 

GC TO ( 1CCC, LCCiflOCl ) , NP R I NT 
100U WRITE (6*11) 

11 FORMAT (//7X4FP/RB, I2X4HZ/RB, 12X4HS/RB*8X1 1HUMEGA THETA* 7X8HCMEGA X 
1 I , 7X10HCE OMEGA X I , 8 X5 HGAMMA/ ) 

Z ED=0 . 

5=0 • 

DO 8 1=1, NMAX 
I E ( 1-1 )8,8,5 
V DEMODE L 

IF ( ( I . EC. 2 ) -OR . ( I .EQ.NMAX ) I CEM = .5*DEL 
5 = S + D E M 

ARGU=UEM**2-(R(I)-R(I-i))**? 

I E ( ARGU .LE .C. ) GC TC 8 
ZED = SURT ( ARGU ) +ZED 

o WR ITF( 6 ,12 )R ( 1 ) ,ZED,S ,CMT ( I ) , OMXI ( I ) , DECMX ( I ) , G AM ( I ) 


12 FURMAT ( 7 E 1 6 • 8 ) 

RETURN 

1001 IF (NPR1M.FC.3) GC TO 51 * 

NMAX1=NMAX-1 * 

DC 6 J= 1 , 4 * 

PEE(J*3,L)=.5*(PFF(J,3,1)+PEE(J,3,2)) * 

PEE { J, 3,NMAX ) = .5* (PEE ( J , i.NMAXi )+PEE( J, 3,NMAX ) ) * 

ZDOT( J ,1) = .5* ( 2 COT ( J, 1 J+ZDCT (J,2) ) * 

ZOO T ( J ,NMAX ) = .5* ( ZDCT ( J , NMA X 1 ) + ZDUT ( J , NMAXi ) * 

ZODUTt J,1 ) = . 5 * [ 7 D CC T ( J , 1 ) + Z 0 DO F ( J , 2 ) ) * 

o ZCOUT ( J , NM A X ) = .5 * ( 7 COU T ( J , N M AX 1 J + ZDOUT ( J, NMAX ) i * 

*KITE(6»6C2) * 

602 FORMAT ( //2EF INITIAL 0 I S P L A C J ME NT S GIVEN// 

14 X l F*5 » 8X4 EL XI,11X7HU T HLT A , 10X 4HM XI,14X1HW/) 

KCUNT=C 

DO 603 1=1, NMAX 


IF ( ( I . EO. 1 ) . CR . ( I.EQ.2 ) .UR. ( I. EQ.NMAX ) ) GO TO 604 
KUUNT=KCUM+ I 

IF (KOUNT-FREC) 603,605,605 
606 KOUNT=C 

604 WR I TE( 6 ,6C6 ) ESS{I) f PEE(l,3,I),PEE(2,3,I),PEE(4 f 3,l),PEE<3,3*l) 
606 FORMAT ( lXFt-3 ,4E16. 8) 

603 CONTINUE 

PRINT 2 CO 1 

2001 FORMAT ( 1H 1 / /* Z CCT*/ / I 
KCUNT=C 

DL 30)0 1=1, NMAX 

IF ( ( I . TQ. I) .CR .( I .FQ.2) .OR. ( I. EQ.NMAX ) ) GO TO 1004 
KCUNT=KCUNT + 1 

IF (KULNT-f PEG ) 3 CCO , 1 0 C6 , I 0 C5 

1006 KCUNT=0 

1004 PRINT 1 C07 , ( l CCT < K , I ) , K= 1 , 4 ) 

1007 FORMA! ( 4E2C .8 ) 

300o CONTINUE 

PRINT 1CC8 


I 00 6 FORMAT (///♦ ZCCOT*//) * 

K 0 U N T = C * 

DU 2000 1 = 1, NMAX * 

IF ( ( I . EQ . 1 ) .CR . ( I . tQ. 2 ) .OR .{ I. EQ.NMAX ) ) GO TO 2004 * 
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KQUNT=KGGN T ♦ 1 * 

IF IKOUNT-FREC) 2 0 CO, 20 05 , 20 05 * 

2005 KCUNT=0 * 

2004 PRINT 1C07 , I Z D DO T (K,I ) , K= I , 4 ) * 

2000 CONTINUE * 

GO TO 6CC 
51 KGUNT=0 

DO 33 I = l,f*PAX 


IF ( ( I.EQ.ll.CP.I I.EQ.NMAX) ) GO TO 103 
14 PEEI4,3,n = XI4,n 
PEEI1,3,I)=XI1,I) 

P EE 1 2, 3 , I ) = X ( 2 , I i 
PEEI3,3,I)=X(3,I) 

GO TO 33 

103 IFII-1) 1 C4 , 1 C4 , 105 

104 K=2 
J=1 

GO TO 1C6 

105 K=NMAX 
J = N PA X 

106 PEE 1 4 , 3 , J )=.5*IXI 4,K)+X(4»K— 1) ) 

PEEt 1 ,3,J )=.5*(X( 1,K)*X (1,K-1)) 
PEEC2,3,J) 5S .5*(X(2,K) + X(2, K— 1 ) ) 

PEE (3,3, J )=.5*(XI 3,K)+X(3,K-1> ) 

33 CONTINUE 


600 JSAVE=JTIPE41 * 

IF IJSAVE.GT.5) JSAVE=5 * 

GO TO 1 1,4,77,10,100, JSAVE * 

1 DO 200 K= 1 , NP A X * 

DO 200 J= 1 , 4 * 

200 ZSAVEI J,1,K) = PEEIJ, 3, K) * 

GO TO 5C1 * 

4 DO 210 K= 1 , NP A X * 

DO 210 J= 1 , 4 * 

210 ZSAVEI J,2,K)=PEE( J,3,K) * 

GO TO 501 * 

77 DO 22 K=1,NPAX * 

DO 22 J=L,4 * 

22 ZSAVEI J,3,K)=FEEI J,3,K) * 

GO TO 501 * 

10 DO 23 K = 1 , A PA X * 

DG 23 J= 1 , 4 * 

23 ZSAVE( J,4,KJ=PEE( J,3,K| * 

GO TO 501 * 

LOO DO 24 K= 1 , NPA X * 

DO 24 J= 1 , 4 * 

ZSAVEI J ,1 ,K) = ZSAVE( J, 2 ,K) * 

ZSAVEI J ,2 ,K ) = ZSAVE I J, 3,K) * 

ZSAVEI J,3,K)=ZSAVE(J*4,K) * 

24 ZSAVEI J,4,K)=PEEl J,3,K) * 

501 IF INPRINT. EC. 2) GO TO 13 * 

IF IMQCtJT IPE,KTHTIPE) .NE. 0) GO TO 1003 * 


NOP TS= INPAX-2 J/FREQ42 
SMAX=DEL*FLCAT INMJX-2 ) 

WRITE! 6,53)SPAX 
53 FORMAT I/9H SP / X/R B= E16 . 8 ) 

E SS 1 1 ) =0. 

DO 20 1=2 , NP AX 
DEM=DE L 

IFI I I. EQ. 2). OP. 1 1 .EC.NMAX) ) DEM=DEL/2. 

20 E SS I I ) = ESSII-D + DEP /SPAX 
601 WRITEI6,2) 

2 FORMAT I/4X1FS,8X4EN XI,11X7HN THETA, 10X4HN XT,12X4HQ XI,12X4HM XI, 
111X7HM THETA/) 

KOU NT = C 

DO 3 1=1, NM AX 

IFI I I.EQ.l ).CR.I I .EC. 2) .OR. I I.EQ.NMAX) ) GO TO 17 
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K0UNT=KClM4l 
I F ( KOUNT-F R EQ ) 3,18,18 
Id KOUNT=C 

17 WRl TEC 6,7 ) ESS(I) , PEE ( 1,2,1 ) ,PEE< 2, 1, I) f PEE < 2 , 2 • I) , P EE ( 3, 2 • I) , 

1 PEE (4* 3, I ) , PEE ( 1 , 1 , 1 ) 

7 FORMAT (lXFt.3 ,6E16.8) 

3 CONTINUE 

MRI TE( 6 »2C 1 ) 

201 FORMAT { /4X1HS,8X4HM XT,12X4HU XI,11X7HU THET A , 1 2X 1HW , 1 2X6HPH I XI/) 

KOUNT-O 

00 202 1 = 1 , KM AX 

IF ((I.EQ.1).CR.II.E0.2).0R.(I.E0«NHAXII GO TO 203 
K0UNT=KCUM + 1 

IF ( KOLNT-F REC ) 202 ,204,204 

204 KGUNT= C 

203 WR I TEC 6 , 20 5) F SS (I) , PEE (3 , 1 , I) , PEE ( 1, 3 , I ) , PEE ( 2 , 3, I ) , PEE ( 3 f 3 , I ) , 

1PEE (4, 2 , I) 

205 FORMAT (1XF6*3,5E16.8) 

202 CONTINUE 
1003 CONTINUE 

500 M=0 

JJ( 1) =1 
J J C 2 ) =2 
JJ( 3) =2 
JJ(A) =3 
J J ( 5 ) =4 

J J ( 6 ) =1 
J J ( 7 ) =3 
J J { 8 ) - 1 

J J < 9) =2 
JJ( 10) =3 
JJ( 1 1 ) =4 
K K ( 1 ) =2 

KK ( 2 ) =1 
KK ( 3 ) =2 
KK ( 4 ) =2 
KK ( 5 ) =3 
KK ( b ) =1 
KK ( 7 ) =1 
K K ( 8 ) =3 

KK ( 9 ) =3 
KK ( 10) =3 
KM 11 ) = 2 
KKK= 1 1 

DC 50 1=1, KKK 
J=J J(L ) 

K= K K ( L ) 

KOUNT=0 

DC 19 I = 1 , K M A X 
M= I 

IF | ( I . t«.l ) .CR.( I .EQ.2 ) ) GO TO 19 

M= ( fSMAX-2 J/FREC+2 

IF ( I «b Q • NM A X ) GC TO 19 

K GU NT-KCUM + 1 * 

IF(K0UNT-FPEC)19,21,21 * 

21 KCUNT=0 

M=( 1 — 2 )/FREC+2 
19 YORDC M ) = PEE C J ,K, I ) 

Nl= ( NM AX-2 J/FPEC + 2 
101. FORMAT (5E12. 5,4X2 14) 


DO 102 1=1, M * 

ABC=YCRO(I) * 

CALL R ECOU T (14, 1, C, ABC) • 

102 CONTINUE * 

50 CONTINUE 
13 RETURN 
END 
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SUBROUTINE EF G t K , l ND5 , NMAX ) 

C SUBROUTINE EFG THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE E , F , AND G 
C MATRICES, AS DEFINED IN APPENDIX A OF REFERENCE ( 1 ) , AT THE STATION SPECIFIED 
C BY THE INDEX K. 

COMMON 

1/BL1/RI 1C2) ,GAMt 102) ,OMTllC2) ,OMXI < 1 02 ) , Dt OMX ( 1 02 ) 

2/BL2/E (4»4),G(4,4)fF(4,4) 

3/BL 3/NU,LAM ,N,EALSIG,CFAR, DEL 
REAL NU , L AM , N , LAM2 . 

CALL BDB(K,CEL,NU,BfDBtD»DD) 

E ( 1 , 1 ) = B 
E<1,2)=0. 

E(1 , 3 ) = C • 

Ed,4) = 0* 

E ( 2 , 1 ) = C • 

D 1= ( 1 * — NU ) 

LAM2- L A M* * 2 
RA^RtK ) 

GA=GAM( K ) 

0X = CMX I (K ) 

0 T=OMT (K) 

DEX = DE C MX ( K ) 

G A2=G A**2 
RE X= ( 3 • *OT“CX ) 

RXE=t3**CX-CT ) 

0TX=0T*0X 

DNLR=L AM2*C*N*C1/ (2.+RA ) 

DONLR=DNLP*CD/D 

E<2,2)=B*Cl/2.+LAM2*0*01*REX**2/8. 

E ( 2 , 3 ) =DNLP*REX 
E(2,4)=C. 

E ( 3 « 1 ) =C • 

E< 3 ,2 ) =E ( 2 ,3) 

K AN= ( N / R A ) **2 

E (3 ,3 ) =LAM2*C*D1* (2.*RAN+( 1 -+NU )*GA2) 

E<3,4)=LAM2 
E ( 4 , 1 ) =C • 

E(4,2)=C. 

E ( 4 , 3 ) =-D 
E(4,4)=C. 

F ( L,1)=GA*B+Ce 

Ft 1 ,2) = ( 1 .+NU)*B*N/ (2 .*RA)+DNLR*REX*KXE/4. 

F ( I ,3 ) = ( CX+NU*OT )+L AM2*D*C 1* t ( 1 • + NU ) *GA2*GX+RAN*RXE/2. ) 

Ft 1 ,4)=LAM2*CX 
F < 2 , 1 ) =— F t 1,2) 

F ( 2 , 2 ) = t Cl /2 . )*(GA*e + OB )-( L AM2* 0*D 1 *R E X / 8 . I * t 2 . *DE X-GA* ( 5 . *CX 
L-3* *0T ) )+LAM2*CD*Dl*REX**2/8* 

F(2 ,3 ) = DNLP*(2 .*( l.+NU ) *GA *0 T-DEX + 3 **GA * ( OX-O T ) ) +DDNL R*REX 
F(2,4)=0. 

F(3,l)=-F< 1,3) 

F(3,2)=DM_R*(3.*GA*UX-GA*0T*( 5 . +2 . *NU ) - DE X ) +DD NLR*R E X 
F < 3 , 3 ) =— LA M2*D*D1* ( (1 « + NU ) * ( 2 . * GA* 0X*0T +GA* * 3 ) + 2 • *GA* R AN ) 
l + LAM2*CD*Cl*( t l.+NU)*GA2+2,*RAN) 

F t 3 ,4 ) =LAM2*GA*(2 .-NU) 

F«4,1)=D*G> ' 

F(4,2)=0. 

F t 4 , 3 ) =— D * NU* G A 
F(4,4)=C. 

Gt L ,1 ) = NU*CB*GA-NU*B*CJTX-B*GA2->0L*B*R AN/2.-L AM2*Q*DI* ( ( 1. + NU ) *GA2* 
lOX**2+RXE**2*RAN/6. ) 

Gt 1 ,2 )=NU*N*0e/RA-(3.-NU)/ < 2 . *R A ) * GA*B * N-DNL R*2 . *G A * ( R EX*R X E / 8 . 

1+t L. + NU)*GT>) 

Gtl ,3 ) = B*t CEX+GA*(CX— OT ) ) + DB* ( 0 X+NU*OT ) -L AM2*D*0 1 *GA*R AN* ( RXE/2. + t 
11* + NU ) *CX ) 

Gt I ,4) = LAM2*D1*GA*0X 

G ( 2 , 1 ) =— B * G A * N * ( 3 • - N U ) / < 2 . * R A ) - D 1 * N *DB / ( 2 . *R A ) +DNLR *2 . * ( - 1 . * t 1- + 

1NU) *GA*0TX+GA/e**t6.*0TX-7.*0X**2-3.*UT**2 )-DEX/4.*{ 5.*0T-3.*CX) ) 
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2-0DNLR/4.*PEX*RXE 

G(2,2)=-C-A*F ( 2,2)+Cl/2 .*B*C T X-B *RA N-L AM2*D*D 1* ( { 1 . + NU ) *0T* *2*R AN 
1— UTX/8 .*RE>**2 i 

G(2 f 3> C7*Nb*CX) /RA+DNLR*(GA*DEX-2.*GA2*OX-2.*( l.+NUI+OT 
1*RAN+REX*{ GA2+CTX) ) — DONLR*R E X*G A 
G ( 2 ,A)=-MU*tAR2*GT*N/RA 

G( 3, 1) =-6*CA* <CT+NU*OX I+LAM2*D*D1* (GA*( U + NU ) * ( -G A*D E X + GA 2*0 X 

1- OX*RAN + 2.*CTX*CX )+FAN/2.*< G A*0 X-G A*GT- 3. *DE X ) ) 

2- LAM2*CC*C 1* ( ( 1 , + NU)*GA2*GX + RAN/2.*‘RXE) 

G ( 3 » 2 ) — — B + N*( CT+NU*CX ) / R A+ DN LR* ( 2 . * ( 1 . + NU ) * ( CT X*OT-G A2 *OX + 2 . *GA2 
1*0T— GT*RAN )+GA*CEX+3.*GA2*( OT-OX) +QTX*REX I -DDNLR*( 2. * { l.+NU)*GA 
2*CT+GA*REX ) 

G ( 3 » 3 ) =“B * ( CX**2+2.*NU*UTX+QT**2)+LAM2*U*DL*RAN*{ { L • + NU ) * ( QTX— RAN 
1+2. *GA2 ) +2 . * ( G A2+OTX ) ) - L AM 2 * CD*D 1 *R AN* ( 3.+NU) *GA 
G ( 3 t 4)=-LAN2*<Dl*CTX+NU*RAN ) 

G(4fl)=C* (C£X3MJ*GA*0X) 

G{4,2) =C*NL*N*CT/RA 
G(4,3) = C«NL*RAfs 
G ( 4 f 4 ) •=— 1 . 

RETURN 

END 


SUBROUTINE FCRCt (K t IN05 ,NMAX) 

C SUuKUUTINE t-URCr THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE LOWER CASE 
C E-VECTJR AS ULFINEC IN APPENDIX A OF REFERENCE ( 1 ), AT THE STATION SPECIFIED 
C BY THE INDEX K . 

C CM RUN 

i/BLl/K ( 1C2) ? G A M ( lC2)tOMT( 1C? ) ,0MXI ( 102 ) ,OEGMX< 102 ) 

3/BL 3/NUtLAR ? N ,EALSIG*CHAR f DEL 
4 /BL 4/C E E ( 4 ) 

REAL NU,LAR,N,L2 
HEAL MSIP 
R A = R ( K ) 

G A = GAM { K J 
JX=CMX l(K) 

U r = CMT ( K ) 

T = T EtMP ( 8 ,CfcU 
DT=DTERP( K f CEL ) 

Of L T 1 =D E L T (K f CEL) 

DLT1=DDELT (K,CELI 
PX1=PX ( K, LFL ) 

P T 1 -P T ( K,CtL ) 

P 1 = P ( K ) 

H = H H T ( K t CE L ) 

U h = DHHT ( K , CtL ) 

HR B = HR A (K , L f L ) 

UHR B= OH R A { K T C EL ) 

Dl= 1 .-NU 
L2=LAM**2 
t AL^tALSIG 

TSUBl = 2.*H*EAL/( D l *HRB ) *T 


MSTP=CHAR* EAL/<3.*C1»*(U .5/HRB-3 . / HRB* *2 + 2 . / HR B**3 >*(DLTl*H**3+3. * 

l *H«*2*CF*UELTl) +DFI. T L *H ** 3 * 0 HR B /HR B ** 2 * ( - 1 . 5 + 6. /HR B- 6 . /HRB * *2 ) ) * 

CtE {4} =CHAR*fcAL*DELTL*H**3/ < 3.*IH ) *( 1 . 5 / HRB- 3 . /HR B**2 + 2 . /HR B** 3 ) * 

CEE( 1)=-PX 1+2.*EAL/(D1*HRB)*(H*CT+T*DH)— DHRB/HRB*TSUBT— * 

1L2*D1*GA*CX*CFE (4 ) * 

CEE ( 2 ) = -PT 1-N/RA*7SU3T-L2*D1*N/RA*0T*CEE(4) * 

CEE I ^)=-Pl-(C X+OT ) *TSUBT-L2*D1*GA*MSTP+ * 

lL2*Ci*Cf E (4 )* <CX*GT-(N/RA>**2) * 

kETURN 
END 
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SUBROUTINE AeCG (K) 

SUBROUTINE AECG — THIS SUBROUTINE CALCULATES THE A,B,C, AND LOWER 
CASE G MATRICES USING THE CURRENT VALUES OF THE E,F,G, AND LCWEK CASE 
E MATRICES. 

COMMON 

2/BL2/E (4,4),G(4,4),F<4, A) 

3/BL3/NU ,LAfr^ ,EALSIG,CFAR» DEL 
4/ BL 4/C EE { 4 ) 

6/BL6/A (4,4),e<4,4) ,C< 4,4) ,SMAG(4) 

9/BL9/D<4,4) ,ZC0T(4,102 ) , ZDDCT ( 4 , 102 ) , EPS 

A/ BL 10/ JT I M E 

B/BL li/ZSAV E (4 ,4,102 ) 

DIMENSION TERMI4I 
REAL N , NU , t A M 
02=2 ./DEL 
04=4. /DEL 
DX=2.*DEL 
DO 11=1,4 

SMAG ( I ) = CX*CEE ( I ) 

DO I J = L ,4 

B ( I , J)=-D4*E ( I, J )+CX*G( I, J ) 

C(I,J>=C2*E(I,J)-F(I,J) 

1 A ( I , J)=C2*E( I , J)+F< I, J ) 

JJ= JTIME+1 

IF (JJ.GT.5) J J=5 
GO TO (2,3,4, 5,6 J , JJ 

3 DO 14 J= I , 4 

L 4 TERM( J)=— 6.*ZSAVE(J ,1 ,K )-6 .* EP $*ZDU T( J , K 1-2. *EP S**2*Z DOO f ( J , K I 
DO 15 1=1,4 
SUM = 0. C 
DO 25 J= L , 4 

25 SLM=SUM + D< I , J)*TEPMJ) 

15 SMAGC l ) = SMAG( l ) + SUM*DX/EPS**2 
DO 16 1=1,4 
DO 16 J= 1 , 4 

lb B( I , J) = ei I , J )-6.*DX*D( I ,J ) / EPS**2 

2 RETURN 

4 DO 17 J = 1 » 4 

1 7 TERM( J)=-4.*ZSAVE ( J ,2 , K )+2 . *ZS AV E ( J , 1 , K )-EP S **2 *Z DDQT ( J , K I 
DO 18 1=1,4 
SUM=0 • 0 
DO 26 J=1 ,4 

2 0 SLM=SUP+D( I ,J)*TERMJ) 

Id SMAG ( I ) = SMAG( I >+SLP*DX/£PS**2 
GO TO 15 

5 DO 20 J=1 ,4 

2 0 TERM(J)=-5.*ZSAVE(J,3,K ) + 4«*ZSAVE(J,2,K)~ZSAVt(J*l,K) 

DO 21 1=1,4 
SUM=0. C 
DO 27 J= 1,4 

27 SLM=SUM+D( I, J)*TERM(J) 

21 SMAG( I ) = S P A C- ( I MSUM*0X/EPS**2 
GO TO 19 

b DO 22 J = 1 » 4 

2 2 T ERM ( J )=— 5 .*ZSAVE(J,4,K )+4.*ZSAVE( J, 3,K |-ZSAVE(J,2,K) 

DO 23 1=1,4 
SUM=0 • C 
DO 28 J = 1 ,4 

2b SLM=SUM + D( I ,J)*TERM(J) 

2 3 SMAG( I J=SMAG( I) +S UM*DX / EPS * *2 
19 DO 24 1 = 1 ,4 
DO 24 J = 1 » 4 

24 B ( I , J) =B( I,J)“2.*CX*D( l , J ) / EP$**2 
RETURN 
END 
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SUBROUTINE I N l T ( INC5 ,BCUNDRY) 

CALCULATION OF TFF P-MATRIX AND THE X-VECTOR AT THE FIRST STATION. 
INTEGER BCLNCPY 
COMMON 

2/BL2/E(A,A) »G (A ,4 ) , F( A , A) 

3/BL3/NL ,LAN,N,EALSIG,CHAR, DEL 
A/BLA/CEE (A ) 

5/BL3/H (A, A > ,FF (A) , JAY ( A ,4 ) 

6/BL6/A(4,4),f!<4,4),C{4,4),SMAG<4) 

7/BL 7/PEE (4,4, l C2 ) * X ( 4 , i 02 » 

DIMENSION CNtGl(4,4),CAPLl(4,4) , EL 1 (4 ) 

U I MENS ICN IPIVCT (A) f INDEX! 4, 2) , A0( A, A) , 80 ( A, A) , A3( A, AI f AA( A, A) t G0( 
IA) 

NAMELISI/F IRST/CMEG1,CAPL1, ELI 
REAL JAY,N,NU,LAM 

90 DO 91 [ = 1,4 

91 GO ( I ) =0 . 

IF ( IN03) 13, 13, 1C 
lu IF (INC5-2) 11,11,13 

11 CALL PCLE(N,DEL,F,G) 

IF (BuLNDRY.EC.il GC TC 16 
*RI Tt< 6,166) 

l6o FORMAT (/27F ECLNCAR Y CCNOITICNS AT S=0 ) 

WR I T E ( 6 , AO ) 

FORMAT (38F CONDITIONS FCR fi SHELL POLE CENERATEO) 

GO TO 16 

1^ IF ( OOUNORY.EC.l ) GO TO TO 
DO Al 1 = 1, A 

eli m =c.c 

DO Al J = 1 » 4 
OMFGK I , J)=C.C 
Al CAPLK I,J)=O.C 
RtA0(3, FIRST) 

13 wKlTEC 6,166) 

WR 1 TE( 6,16?) 

167 FORMAT ( /2 l X C M FGA , 49 X 5HL A MD A , 3 4X 3 HEL L / ) 

DO 16 8 1=1, A 

16 6 wk I Tt( 6 ,16S KMEGl < 1,1) ,CMEGl ( I , 2 ) , G ME G 1 (I , 3 ) , OM EG l < I , A ) , 

1 CAPL1 ( I ,1 ) , CAP L 1 II, 2) ,CAPL 1 ( 1 , 3) , CAPLK I , A) , 

2 EL 1 ( I ) 

Lo9 FuRMAT (4c 12.4,2(6X4612.4) ) 

70 DO l 1=1, A 
DO 1 J=l,4 

B0( I ,J)=h( I ,J)/CEL*JAY( I , J ) / 2 • 
i AU(I,J)=JAYl I , J ) / 2 . -H ( I , J ) / C t L 
DC 2 1=1, A 
DO 2 J = 1 , A 
S l = 0 . 

S2 = U . 

DO 3 L= 1 , A 

Sl=$i+CNtGl (I ,LJ*AO(L, J) 
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3 S2=S2+GPEG1 (I ,L)*80(L* J> 

A3 ( I f J ) = S 1 

2 A4( It J)=S2 
DO 5 1=1*4 
S1 = 0. 

00 4 J = 1 * 4 

S1=S1 + 0MEG1< I * J)*FF(J) 

BO ( 1 1 J )=A3 ( I*J)+CAPL1( I * J ) / 2 • 

4 A0( I, J)=A4I I* JI+CAPLII 

5 G0( I ) = ELI ( 1 )-S 1 

CALL MATINV ( C ,4 , CE E, 0 * DETE RM , I PIVOT* I NCEX*4* I SCALE ) 
DO 50 1=1,4 
DO 50 J = 1 * 4 
S 1=0. 

S2= 0 • 

S3= 0. 

00 51 K = 1 , 4 
S1=S1+C(I,K)*E(K,J) 

S2= S2 + C ( I ,K)*A(K,J) 

51 S3=S3+BC( I,K)*CfK,J| 

A 3( 1, J) = S1 

A4{ 1, J ) =S 2 
50 E(I , J)=S3 
DO 52 1=1*4 
DO 52 J = 1 * 4 
S 1= 0 . 

S2= 0 • 

DO 53 K = 1 ,4 
S 1= S 1 + B 0 ( I ,K)*A3( K, J) 

53 S2=S2+B0( I »K)*A4( K*J) 

F< I*J) = S1“AC(I *J) 

52 G ( I , J)=S2 

16 CALL MA7INV(F,4,G,4,DETERM, IPIVOT, INDbX,4, ISCALt) 

IF( IN05.EG.C.CR.IND5.EQ.3) GO TC 59 
DO 60 1=1,4 
Xt I * 1 ) = C • 

DC 60 J = 1 * 4 

60 P EE ( I * J , l ) = G ( I*JJ 
CALL P A ND X ( 2 ) 

RETURN 

59 DO 54 1=1,4 
S1=0. 

DO 61 J = 1 ,4 
P EE { I * J ,2 ) = G ( I,J) 

61 S1=E( I,J|*SPAG(J)+S1 

54 GO ( I > = S1-GC { I ) 

DO 56 1=1,4 

S 1 = 0 • 

DC 57 J = 1 ,4 
57 S 1= S 1 + F ( I , J)*GC( J) 

56 X< I ,2 ) = S 1 
RETURN 
END 



APPENDIX D - Continued 


SUBROUTINE PANCX(K) 

C SUBROUTINE PANCX THIS SUBROUTINE CALCULATES THE P MATRIX AND THE X-VECTOR 
C AT THE STATION SPECIFIED BY THE INDEX K USING EOUAT ION (29) OF THE TEXT. 
CCMMON 

6/BL6/A (A* A > ,B(A,4) ,C< 4 , 4 ) , SMAGt 4 ) 

7/BL7/PEE<4,4, 1C2),X(4, IC2) 

D I MENS ICNP 1 ( A * A ) » I PIVOT (A) , INDEX (A, 2) t X 1 ( A ) » X2 ( A ) 

DO I 1 = 1,4 
DO 1 J = 1 , 4 
SUM=0 • 

DC 2 L= 1 * 4 

2 SLM=SUM*C( I,L)*PEE(L,J»K— 1) 

L Pit I* J ) = B( I , J ) — SUM 

C ALL MAT IN V (Pi , 4 , X2 , G , D E TER M , I P I VOT , I ND E X , 4 , I SC AL E ) 

DO A 1=1, A 
SUM=0 • 

DO 3 J = 1 , A 

3 SUM=SUM+C(I,J)*X( J,K-1) 

A Xl( I ) = S MAG ( I J-SLM 

DO b 1=1, A 
DC 5 J = 1 , A 
SUM=0. 

DC 6 L= 1 i A 

o SUM=SUM+P1 (I * L ) * A ( L , J ) 
b P t E ( I , J , K ) = SU M 
DO 7 1=1, A 
SUM=0 • 

DO 8 J = 1 , A 

b S LM = SUM + P 1 ( I » J ) * X 1 ( J) 

7 X( I ,K)=SUM 
RETURN 
END 


SUBROUTINE EC73(K) 

C SUBROUTINE EQ73 THIS SUBROUTINE CALCULATES THE SOLUTION VECTOR AT 
C THE STATIUMK) , GIVEN THE SOLUTION AT K+ 1 . 

COMMON 

f/BL7/PEE(A,A,lC2) ,X(A, 102 ) 

DIMENSION iMA,102) 

EQUIVALENCE ( X ( l , 1 ) , Z < 1 , 1) ) 

DU l 1 = 1, A 
S UM = U • 

DU ? J= 1 , A 

2 SUM=SUM + PEf { I, J,K )*Z( J , K+ 1 ) 

1 Z(I,K) = X(I,K ) - SLM 
RETURN 
END 
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SUBROUTINE FINAL ( N MAX • IND5 , BOUNDR Y ) 

C CALCULATION CF SOLUTION VECTOR ASSOCIATED WITH LAST STATION 

INTEGER BCUNCRY 
COMMON 

3/BL3/NU,LAM, N, EALSIG,CHAR, DEL 
5/BL5/H (4,4 ) ,FF (4) , JAY(4,4) 

7/0L7/PEEI 4, 4, 1 C2 ) « X ( 4 , 1 C2 » 

DIMENSION CMEGL(4,4),CAPLL(4f4) *ELL(4) 

DIMENSION I PI VCT ( 4 1 , 1 NO EX ( 4 f 2 ) , A1 ( 4, 4) , A2 (4 , 4 ) * A3 (4 f 4 ) , PS I ( 4 , 4 ) f 
lGM(4 f 4) f ETA (4),B(4,4) 

REAL JAY , N , NU , L AM 
NAMELI ST/LAST/CMEGL,CAPLL, ELL 
90 IF C IN 05 I 13,13,10 

10 IF ( IND5-2) 13,11 ,11 

11 CALL P CL E IN, CEL, PS I,GM I 
IF (BULNCFY.EC.l) GC TO 16 
WRITE(6,17C) 

170 FORMAT ( /30F BCLNDARY CONDITIONS AT $=SMAX> 

WR I T E ( 6 ,4 C ) 

40 FORMAT (38F CCNCITIONS FOR A SHELL POLE GENERATED) 

GO TO 16 

13 IF IBOLNDRY.EC.l) GO TO 70 
DO 41 1=1,4 

E LL ( I ) =0 • 0 
DO 41 J= 1 , 4 
OMEGL ( I , J ) = 0.C 

41 C A P LL ( i , J ) = C.C 
READ IS, LAST) 

13 WRITE! 6 , 1 7 C ) 

WR I T El 6 ,167 ) 

16 J FURMAT !21XSFCMEGA,49X5HLAMDA,34X3HELL/ ) 

DO 171 1=1, A 

171 WRITE! 6 , 1 6 S ) CMEGL (1,1) , CMEGL ( I , 2 ) , OMb GL (I , 3 ) , OMEGL! 1,4), 

1 CAPLL ( I ,1 ) , CAP L L ( I , 2 ) , C AP L L < I , 3 ) , C A PLL ( I , 4 I , 

2 ELL(I) 

169 FORMAT (4E1^*4»2(6X4E12.4) ) 

70 DO 1 1=1,4 
DU 1 J= 1,4 

Alt I , J > = JA Y ( I,J)/2*+H( I fJ)/DEL 

1 A2! I, J) = JAY(t,J)/2.-HU,J)/DEL 
DO 2 1=1,4 

DU 2 J = l , 4 
S2=C. 

S3= 0 • 

DO 3 L= 1 ? 4 

S2= OMEGL ( I,L)*A1(L,J)+S2 
3 S3 = UMEGL! I ,L)*A2(L,J)+S3 
PSI (I , J)=S3+CAPLL ( I, J) /2. 

2 GM( I,J) = S2+CAPLL( I , J ) / 2 • 

16 DC 4 1=1,4 

DC 4 J= 1 , A 
S1 = 0. 

DU 5 L = 1 , 4 

3 Sl = Si + PSI( I,L)*PEE(L,J,NMAX-1) 

4 B ( I , J ) = GM ( I,J)-S1 
DO 6 1=1,4 

S1 = 0, 

S2=0. 

DO 7 J= 1,4 

S1=S1*PSI ( I , J ) *X ( J , NMAX-I ) 

7 S2=S2+GMEGL ( I , J ) *FF (J) 

6 ETA! I ) = b L L ( I J-S1-S2 

CALL MATINV (E ,4* ETA, 1, DETERM, IPIVQT , INDEX, 4, I SC ALE) 

DO 12 1=1,4 
12 X(I , NM A X) = tTA(I) 

RETURN 

END 
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SUB ROUT INt $TRFSS(FREQ,NMAX, I ND5 ) 

SUBROUTINE STRESS — THIS SUBROUTINE CALCULATES THE SECONDARY QUANTITIES 
N XI, N XI THETA, Q XI, PHI, M THETA, M XI THETA, M XI THETA, AND 
N THETA AT EACH STATION ALONG THE SHELL • 

COMMON 

I/BL 1/R ( 1C2 ) , G AM ( 1 C2 ),GMT ( 102 ) »OMXI ( L02 ) ,DECMX< 102) 

3/BL3/NU ,LAH,N,EALSIC,CHAR, CEL 
5/BLb/H (4,4) , F F ( 4 ) , JAYl 4,4) 

7/8L7/PEE(4 ,4, IC2 ) ,X (4, 1 C2 ) 
d/BL8/AK(3 ,4 ), ALL (3,4) ,STHER(3> 

DIMENS ICN l (4, IC2 ) , Y(4 ) ,DZ (4 ) 

EQUIVALENCE (Z(I,1),X(1,I)) 

INTEGER FREC 
REAL N,NU,JAY,LAM 
KCUNT=0 
DU 9 I = 1 , N PA X 

IF( < I .EC. 1 ) .CP. (I -EC.NMAX) > CO TO 1 
I F ( 1-2) 13,13,14 

14 KUUNT=KCUM + 1 

IE( KGUNT— FREQ 19,12, 12 

12 KUUNT=C 

13 DC 3 L= 1 , 4 

Y ( L ) — Z ( L , I ) 

3 DZ ( L )= ( Z (L , 1+1 l-Z IL, I- L ) J/2./DEL 

GO TO 2 

1 IH 1-1 )4,4 ,5 

4 K=2 

GO TO 6 
b K=N M AX 

6 DO 11 L = 1 ,4 

Y { L ) = • 5 * ( Z l L , K ) + Z(L ,K-1 ) ) 

11 DZ ( L) = (Z( L ,K >-Z ( L ,K-1 )) /DEL 

2 CALL HFJ( I, INC5,NMAX, 1. > 

CALL KLT ( I , INCb , NR AX ) 

DO 7 L= 1 ,4 
S U M 1 = 0 . 

S U M 2 = 0 . 

008 M= 1 , 4 

SUM1=SUR1+ H(L,R)*DZ(M> 
d SUM2 = SLR2 + JAY( L , M ) * Y ( M ) 

7 PEE (L, 2 , l ) = SUR1+SUM2 + FF (L) 

DO 9 L = 1 , 3 

SUM 3 = 0 . 

SUM 4= C • 

DO 10 M = 1 , 4 

SLM3 = SLM3 + AK(L,M*CZ(M) 

10 SUM4=SUM4+ALL (L,M)*Y(M) 

PtE (L, 1 , I )=SUR3+SUM4+STFER (L) 

9 CONTINUE 
KETURN 
END 
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SUBROUTINE KLT ( K , I ND5 , NMAX ) 

C SUBROUTINE KLT THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE MATRICES 
C WHICH ALLOW THE CALCULATION OF THE QUANTITIES M-THE TA ,N- THETA , AND M-XI THETA 
C AT THE STATION SPECIFIED BY THE INDEX K. 

COMMON 

I/BL I/R < 102 ) ,GAM( IC2 ), OMT< 102 ) , OMXI ( 102) , DEQMXC 102) 

3/BL3/NU,LAM,N,EALSIG, CHAR, DEL 
8/BL8/AK (3 *4 ) , ALL ( 3 ,4) « ST HER ( 3) 

REAL NU » L A ► » N 

CALL BOB! K,DEL»NU,B,DR,D,DD) 

D1^D*< 1 .-NU ) 

D2=D* ( 1 NL**2 ) 

GX=QMX I (K ) 

eal=ealsig 

H-HHT ( K , D E L ) 

HRB = HR A ( K , C EL ) 

TEMPER=TEMF(K,CEL) 

DELTAT = CELT ( K , CEL ) 

STHER ( 1 )=-CFAR*EAL*CELTAT *H**3 * ( 1 . / ( 2 . *HRB >-l./HRD**2 
l+2./(3.*HRE**3) ) 

STHER ( 2 )=-2.*F*EAL*T6MPER / i ( l.-NU) *HR B ) 

STHERI 3 ) = G • 

I F ( IND5 )1, 1,2 

2 IF I < < I NC5-2 ) .LE.O ) .AND. (K-EC .1 ) ) GO TO 8 
IF l I ( IND5-2 ) .GT.O ) .AND. IK. EQ .NMAX ) ) GO TO 8 
1 R A=R { K ) 

GA= GAM ( K ) 

0T= GMT ( K ) 

REX=3.*CT-CX 
RXE= 3 . *CX— C T 
RAN=N/RA 
R AN2=R A N* * 2 
AK ( 1 v 1 )=0. 

AK ( 1,2 ) =0 . 

AK ( 1, 3)--GA*D2 
AKC 1,41=0. 

AM 2, 1 )=B*KL 
AK( 2,2 ) = 0 . 

AK (2,3 ) =C . 

AK (2,4) =0 . 

AK( 3,1 ) =0 . 

AM 3,2 ) = D 1 *R E X /4 . 

AKI 3,3 ) =RAMCl 
AK ( 3, 4 ) = C • 

ALL ( 1 , 1 ) = G A *C X + C2 

ALL ( 1 * 2 ) — D2*R A K*0 T 

ALL (1,3 >=C2*RAK2 

ALL ( 1 , 4 ) =NL 

ALL ( 2 , 1 ) = B*GA 

ALL ( 2 , 2 ) = 8*R A N 

ALL (2, 3 ) = B*(CT + NU + CX) 

ALL (2,4 ) = C. 

ALL(3, 1 )=-Dl*RAK*RXE/4. 

ALL ( 3 , 2 ) G A* C l + REX/4. ■ 

ALL<3,3)=-GA*RAN*C1 
ALU3,4) = C. 

GO TO 6 
B DO 3 1=1,3 
DO 3 J= 1 , 4 
A K ( I,J)=0. 

3 ALL ( I , J )=C . 

Cl= (N**2/2 .-1 . )/( l.+NU-N**2*NU/2. ) 

AM 1,1 ) =D2*CX* < (l.+NU) *C1 + 1 . ) 

A K { 1,2 )=D2*CX*N*(NU+Cl+1* ) 

ALL (1,4 ) = < AU- <1.-NU**2 >*C1 ) 

AK ( 2, 1 )=B*( l. + NU) 

AK ( 2,2 ) = B*K 
ALL(2,3) = e*(l.-*NUMCX 
C2=2. + 2.*M-NL*N**2 
AK(3,lJ=Din*l 1./C2-OX/2. ) 

ALL ( 3,4)=— N*( 1. — NU) /C2 
6 RETURN 
END 
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SUBROUTINE BDB<K,DEL,NU,B,OB,D,DD) 

SUBROUTINE BCD — THIS SUBROUTINE CALCULATES THE BENDING STIFFNESS 
t Of THE MEMBRANE ST IF FN ESS » B * AND THE DERIVATIVES OF D AND B,DD AND 
DB, RESPECTIVELY, FOR A SHELL COMPOSED CF A CORE HAVING NO STIFFNESS 
AND TWO SYMMETRICAL COVER PLATES 
REAL NLtN»LAM 
HR B = HR A ( K , C EL ) 

UHK B = DHP A ( K,C£L) 

H = HHT ( K , D E L ) 

DH= DHH T ( K , C EL I 

D2=1.-NL**2 

B = 2 • *H/ ( C2 *HR B ) 

D = H **3*< 3 ./ (2 .*HRB )-3 . /HRB**2+2 ./HR 8**3 )/(D2*3. ) 

CB=2.*CE/ ( C2*ERB)- B*CHR B/HRB 

DD=3.*CH*U/H4H**3*DHRB/ (D2 *HRB**2 ) * (- . 5+2 . /HRB- 2. /HR B**2 I 
DD=OU/3. 

RETURN 

END 


SUBROUTINE PCLE(N,DEL,A1,A2) 

C SUBROUTINE POLE- THIS SUBROUTINE CALCULATES THE FINITENESS CONDITIONS FOR 
C A CLOSED SHELL 

OIMENSICN A i ( 4 , 4 ) , A 2 ( A , A ) 

REAL N 
DO I 1 = 1, A 
DO I J= I , A 
A 1 ( I, J)=0. 

1 A2 { 1 , J ) =C . 

I H ( N.EC .0 . > GC TO 2 

IF{ (N.EC. 1. ) .CP. (N.£sJ.-l. ) ) GO TO 3 

AI( 1, 1 ) = . 5 

A 1 ( 2,2 > = .5 

A2< 1,1 ) = .5 

A2 ( 2,2 ) = ■ 5 

IMN.NE.2.) GC TO A 

A2 ( 3 , 3 ) = 1 • / CE L 

All 3,3)=— I ./DEL 

A1 (A, A ) =— 1 . /C E L 

A2( A,A) = I ./CEL 

RE TURN 

A A 1 (A, A ) = • 5 
A2(A,A>=.5 
Ai ( 3,3 )=.5 
A2C 3,3 ) = .5 
RETURN 

2 All 1,1 ) = .S 
All 2,2 > = .5 

A 1 ( 3 , 3 ) ■=— 1 ./DEL 
A 1 ( A , A ) =— 1 • /DEL 
A 2 ( 1, 1 )=.5 
A2( 2,2 ) = . 5 
A 2( 3,3) = 1./DEL 
A2 (A,A)=1./C£L 
. RETURN 

3 All 2, i ) = .* 

A 1 ( 2 , 2 ) = . 5 

All L , 1 ) =- 1 • /DEL 
A L ( 3,3)=. 5 
A 1 (A, A) =.5 
A 2 ( 2, L )=.5 
A2 (2,2 ) = .S 
A2 ( 1, 1 )=1 ./DEL 
A 2 ( 3,3 ) = .S 
A2( A, A )= .5 
RETURN 
END 
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SUBROUTINE INPUT ( NMAX , IND6 ) 

COMMON 

1/BL l/R(102) f GAMlC2),0MT( 1C2),0MXI ( 10 2) ,DE0MX( 1C2) 

3/BL3/NU,LAP,N, EALSIG,CHAR, DEL 

9/BL9/D(4,4) ,ZCGT( 4, 10 2) ,ZDD0T(4,102),EPS * 

6/BL 14/ AO , T I M , C P S * 

REAL N , NU ? LAP * 

PI=3. 14159265358979 * 

RI= 40.1 
R0=90 • 

ANG= PI/6. * 

DEL = (RO- PI) / (RO*CCS( ANG) * E LOA T ( NM AX-2 ) ) * 

R ( NMAX ) = 1.0 * 

R ( 1 ) = RI/RC * 

DELR= ( R ( N FAX ) — R ( 1 ) )/ FLO AT ( N MAX-2 ) * 

R ( 2 ) = R ( 1 ) + DELP/2 * 

N M 1 = NPAX-1 * 

DO 1 1 = 3, NP1 * 

1 R ( I ) = K ( I— 1 i + C EL R * 

DO 2 1=1, PPAX * 

GAMt I )= COS (ANG )/ R ( I ) * 

0 M T ( I ) = SIMANG)/P(I) * 

OMX U I )= 0. * 

2 DEOP X ( I ) = C - * 

C IF RUNTYPE = 2 OR 3 THEN SET EE AND RHO TU FIND TIME INCREMENT, EPS. HERE * 
C EE IS THE PCCULIS CF ELASTICITY AND RHO IS THE DENSITY(L8S/IN**3) * 

GACC= 386.CE8527 * 

t£= I05C0CCC. * 

RHO A= C.l * 

RHQ= RHCA *2.75/27. * 

SS= EE/RFC * 

EORHO= 5 S*GACC * 

LPS= SCPT ( ECRFC/CFAR**2 )*T IM * 

RETURN 
END 


FUNCTION HHT { K ,DEL ) 
HHT = I . 

RETURN 

END 


FUNCTION TEPP(KfDEL) 
T EM P=0 . 

RETURN 

END 
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FUNCT1CN CT EF F ( K t CE L ) 
l> Tt MP = C • 

Kb TURN 
fcfNJ 


FUNCTIUN CbLT{K t OEL) 
D t L T ~ 0 . 

RETURN 
t i\ll) 


FUNCTION CFFT{K t CEL) 
L)HH T =0 . 

RETURN 

END 


FUMCTICN FRA ( K r C E l ) 
H K A = 2 7 * 
l<ET URN 
LNO 


FUNCTIUN CFRA(K,DEL) 
DHRA-0 • 

Rt TURN 
UNO 


E UNCT ICN CCELT t K , CEL) 
ODE L T = C . 

RETURN 

END 
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. FUNCTION P (K) 

COMMON /BL3/NL', LAM, N, EALSIG, CHAR, DEL 

COMMON /BL9/ C(4,4), ZDUT(4,102>, ZDOOT C 4 , 1 0 2 ) , EPS * 

C CM MON /B L1C/ JTIME * 

C/BL12/INITZ * 

6/BL 14/AO, T IM,CPS * 

COMMON /BL16/TAU * 

K/BL18/R0NTYPE , JUMP C 5 ) * 

INTEGER RLNTYFE * 

HEAL N * 

REAL LAM * 

DATA PI/3.14159265358979/ * 

C PS= L 0 * * 

TIM=. 05/CPS * 

TAU= FLOAT (JT IME) *TIM * 

PSIBAR a 5.+PI/180. * 

F EE = P I /6 . * 

PS 1= PSI8AR*SIMCPS*TAU*2.*PI ) * 

Q=l • /LA M * 

SA=SIN(PSI) * 

CA=COS(PSl> * 

SF= SIN (FEE) * 

CF~ COS (FEE) * 

PG= -Q* (SA*SA*SF*SF+ 2 . *CA*C A*C F*CF J * 

Pl= -4.*Q*CA*SA*SF*CF * 

P2— -Q*SA*SA*SF*SF * 

IF (N .EG. C.) P=PO * 

IF (N .EC. 1.) P=P1 * 

IF (N .EC. 2. ) F=P2 * 

RETURN 
END 


FUNCTION P T ( K , CE L ) 
P T -0 . 

RETURN 

END 


FUNCTION PX(K,CEL) 
P X= 0 . 

RETURN 

END 


FUNCTION ZIMTCKK, I) * 

KK= 1 IS FOR ZIM T= U X|(|,0) * 

KK= 2 IS FCR Z I N I T= U T HET A ( 1,0) * 

KK= 3 IS FOR Z IN 1 1= W(I,0) * 

Z I N 1 T= 0 . * 

RETURN * 

END * 
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FUNCTION ZCIN1TIKK, I ) 

KK= I IS FCR Z L IN I T = U XI D0T(I,0) 

KK= 2 IS FOR ZDINIT= U THETA DOT(ItO) 

KK= 3 IS FCR ZCIN IT= M DOT! 1*01 

Z D I M T= 0. 

RETURN 

END 


PKUGKAM SUMUF ( INPUT * OUTPUT, TAPE 1 4 I 
DIMENSION A N S ( 1 1 * 39, 21) 

INTEGER FREC 
INTEGER RUNTYPE 

11 FORMAT ( / 1 > * S 1 A* 2X*N XI*8X*N THETA*5X*N X I THET A*3X* SFEAR*7X 

I *M XI*8X*M TFETA*5X*M XITHETA*3X*U XI *6X* UTH ETA* 6X* W* 1 1 X*PH I * ) 

12 FORMAT ( IX, 12, L1E12.3) 

14 FORMAT <*l*bX*IFE TIMESTEP IS*16,5X*T I ME = *E 15.7) 

UAT A P 1/3. 14159265358979/ 

REWIND 14 

CALL REC1M14, l, IC, FREQ, NMAX , NT t ME , NFOUR , AK, TIM, RUNTYPE) 

L L= 1 

1 F ( RUN TYPE.FC* 2) LL = 2 
THE T-A K *P I 

PRINT b, FREC, NMAX, NT I ME , NFCUR , THET 
:> F CRM AT (bX* FREC = * 16, 10X*NMAX = *I6 , 10 X*NT IME=*I6, 1GX*NF0UR=*I6 
1 /bX*THFTA = < : E 14.7 ) 

KK =11 

N 1 = ( NM AX — 2 ) / F R F Q + 2 
NFULR=NFGLR+1 
DO 20 J=l,Nl 
DC 20 I = L , KK 
DO 20 K = 1,M[ME 
20 ANS(1 , J,K)=0. 

DO 3 M=l, NFCUR 
A N = M— 1 
ANG = AN* THET 
DO 3 K = L L , N T IM£ 

DU 3 1= 1 ,K K 

DO 3 J=1,M 

CALL KtCIM 14, 1, IC, ABC) 

IF ( ( I .Eg* 3) .CR.U.EQ.7) -OR. ( I . fcQ . 9 ) -OR - (I . EQ . 1 1 ) ) GO TO 2 

AN S(I,J,K)=ANS( I , J,K) +ABC*COSt ANG) 

GC TO 3 

2 ANSI I , J»K ) = ANS( I ,J»K) ♦ A PC* S I N ( ANG ) 

3 CLNTINUE 

DO lJ K = LL ,NT IME 
NCYC= K-L 

T 1 M E = T I M*F L C A T INCYC) 

PRINT 14, NCYC, TIME 
PRINT 11 
DO 13 J = 1 , M 
L = J + ( J-2 ) * (FRfcC-1 ) 

IF { J . E Q . 1 ) L = 1 

IF (J -EC- 2) L = 2 

IF (J . EQ.KMAX) L = NNAX 

PRINT 12, L, (ANSI I,J,K), I = 1 ,KK ) 

13 CONTINUE 
STOP 

END UF SUMLP 
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